{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1bf0acd8-c2f1-4ad9-b617-5bdc4262544c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import glob\n",
    "import numpy as np\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import torchvision\n",
    "import math\n",
    "\n",
    "from unet import UNetModel\n",
    "from utils import *\n",
    "from diffusion import GaussianDiffusion\n",
    "\n",
    "import tqdm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e33c8520-d4a0-4ad8-b2f9-61eb7181b222",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e9bcb9ec-fd2b-43c7-9cc8-1d09ba0c8042",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "from torch.utils.data.dataset import Dataset\n",
    "\n",
    "class PatchDataset(Dataset):\n",
    "    def __init__(self, naip_paths_list, label_paths_list, patch_size, mapping=None, transforms=None):\n",
    "        '''\n",
    "        Args:\n",
    "            naip_paths_list (list): Paths to NAIP images\n",
    "            label_paths_list (list): Paths to labels\n",
    "            patch_size (tuple): Extracted patch size\n",
    "            mapping (numpy.array): Label mapping array\n",
    "            transforms (function): Function that applies image and label transforms\n",
    "        '''\n",
    "        self.naip_paths_list = naip_paths_list\n",
    "        self.label_paths_list = label_paths_list\n",
    "        self.patch_size = patch_size\n",
    "        self.mapping = mapping\n",
    "        if transforms is None:\n",
    "            self.transforms = lambda x: x\n",
    "        else:\n",
    "            self.transforms = transforms\n",
    "\n",
    "        # Load tiles\n",
    "        self.naip = []\n",
    "        for path in self.naip_paths_list:\n",
    "            with rasterio.open(path) as f:\n",
    "                self.naip.append(f.read())\n",
    "        print('Loaded', len(self.naip), 'NAIP tiles')\n",
    "\n",
    "        self.labels = []\n",
    "        for path in self.label_paths_list:\n",
    "            with rasterio.open(path) as f:\n",
    "                # Apply label mapping\n",
    "                label_raster = f.read()\n",
    "                label_raster = self.mapping[label_raster.astype(np.uint8)]\n",
    "                self.labels.append(label_raster)\n",
    "        print('Loaded', len(self.labels), 'label tiles')\n",
    "\n",
    "    def __len__(self):\n",
    "        assert len(self.naip_paths_list) == len(self.label_paths_list), 'NAIP/Labels tiles mismatch'\n",
    "        return len(self.naip_paths_list)\n",
    "\n",
    "    def __getitem__(self, idx):\n",
    "        '''\n",
    "        Args:\n",
    "            idx (int): Index of tile to sample from\n",
    "        '''\n",
    "        # Select the tile to get a random patch from\n",
    "        naip = self.naip[idx]\n",
    "        labels = self.labels[idx]\n",
    "        \n",
    "        # Randomly select patch in tile\n",
    "        height, width = naip.shape[1], naip.shape[2]\n",
    "        row = torch.randint(0, height - self.patch_size[0], size=(1,)).item()\n",
    "        col = torch.randint(0, width - self.patch_size[1], size=(1,)).item()\n",
    "\n",
    "        # Use copy to not modify the original labels\n",
    "        naip_patch = naip[:, row:row+self.patch_size[0], col:col+self.patch_size[1]].copy()\n",
    "        labels_patch = labels[:, row:row+self.patch_size[0], col:col+self.patch_size[1]].copy()\n",
    "        \n",
    "        # Apply transforms\n",
    "        naip_patch = torch.from_numpy(naip_patch)\n",
    "        labels_patch = torch.from_numpy(labels_patch)\n",
    "        naip_patch, labels_patch = self.transforms(naip_patch, labels_patch)\n",
    "        \n",
    "        return {'naip': naip_patch, 'labels': labels_patch}\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2b5930f6-19ca-4505-be9c-f2d947648730",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded 10 NAIP tiles\n",
      "Loaded 10 label tiles\n",
      "Created dataset\n"
     ]
    }
   ],
   "source": [
    "# Specify tiles\n",
    "naip_paths = glob.glob('data/enviroatlas_lotp/pittsburgh_pa-2010_1m-train_tiles-debuffered/*naip.tif')\n",
    "naip_paths.sort()\n",
    "label_paths = glob.glob('data/enviroatlas_lotp/pittsburgh_pa-2010_1m-train_tiles-debuffered/*highres_labels*.tif')\n",
    "label_paths.sort()\n",
    "\n",
    "# Transforms to apply\n",
    "def transforms(img, labels, p_horizontal=0.5, p_vertical=0.5):\n",
    "    # Horizontal flip\n",
    "    p = np.random.rand()\n",
    "    if p < p_horizontal:\n",
    "        img = torchvision.transforms.functional.hflip(img)\n",
    "        labels = torchvision.transforms.functional.hflip(labels)\n",
    "        \n",
    "    # Vertical flip\n",
    "    p = np.random.rand()\n",
    "    if p < p_vertical:\n",
    "        img = torchvision.transforms.functional.vflip(img)\n",
    "        labels = torchvision.transforms.functional.vflip(labels)\n",
    "        \n",
    "    # Resizing\n",
    "    img = F.interpolate(img.unsqueeze(0).float(), scale_factor=1/4, mode='nearest').squeeze(0) / 255.\n",
    "    labels = F.interpolate(labels.unsqueeze(0).float(), scale_factor=1/4, mode='nearest').squeeze(0)\n",
    "\n",
    "    return img, labels\n",
    "    \n",
    "# Create data set\n",
    "device = torch.device('cuda:0')\n",
    "patch_size = (256,256)\n",
    "n_labels = len(enviroatlas_simplified_labels)\n",
    "dataset = PatchDataset(naip_paths, label_paths,\n",
    "                       patch_size=patch_size, \n",
    "                       mapping=enviroatlas_simplified_label_mapping_array,\n",
    "                       transforms=transforms)\n",
    "\n",
    "# Create random patch-sampling dataloader\n",
    "batch_size = 32\n",
    "sampler = torch.utils.data.RandomSampler(dataset, replacement=True, num_samples=batch_size)\n",
    "dataloader = torch.utils.data.DataLoader(dataset, batch_size=batch_size, num_workers=0, sampler=sampler)\n",
    "print('Created dataset')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f598c0b4-9aea-4db9-8183-258555bc88a1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAEtCAYAAADHtl7HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABbIElEQVR4nO29d5xdV3X3vdbtc2fu9KJR782yii1LcsG4YzokQCAUJxgcniQEnvA8QJ68pOd9eJOQQBqEYAKBUIwJ2BgwNsLGuFu2Zav3NtL0mTv19rvfP+Zad35LsmaOpl1Jv+/no8/MmnPOPuvss+/WnnN+81vqnBNCCCGEEDJ+fDOdACGEEELIhQYXUIQQQgghHuECihBCCCHEI1xAEUIIIYR4hAsoQgghhBCPcAFFCCGEEOIRLqAIIYSQAqr6qKp+aLqPJRceXECRc6KqR1X1lpnOgxBCvML5i0wlXEARQgghhHiECygyLlT1t1T1CVX9B1WNq+phVb2m8PMTqtqhqneM2v+NqvqiqvYXtv+Zae8DqnpMVbtV9TOjf1NUVZ+qflpVDxW236OqtdN8yYSQixBVrVHVB1S1U1V7C9/PNbstUdVnVbVPVe8bPf+o6hZVfbIwD76kqje8ynmWquovC210qep3p/CyyAzABRTxwmYReVlE6kTkWyLyHRG5SkSWisj7ROSfVbWisO+QiHxARKpF5I0i8j9U9W0iIqq6WkT+VUTeKyLNIlIlInNGnecPRORtIvJaEZktIr0i8i9TdlWEkEsJn4j8h4gsEJH5IpIQkX82+3xARD4oI/NPVkT+UUREVeeIyI9F5K9EpFZE/peIfF9VG85ynr8UkYdEpEZE5orIP032hZCZhQso4oUjzrn/cM7lROS7IjJPRP7COZdyzj0kImkZWUyJc+5R59wO51zeOfeyiHxbRhZEIiLvEJEfOeced86lReRPRGR0UcbfEZE/ds61OOdSIvJnIvIOVQ1Mx0USQi5enHPdzrnvO+eGnXMDIvLXUpybXuEbzrmdzrkhEfmMiLxLVf0y8oviT5xzPynMbQ+LyDYRecNZTpWRkUXabOdc0jn3+NRdFZkJuIAiXmgf9X1CRMQ5Z39WISKiqptV9ZHCY/I+EfmIiNQX9pstIideOcg5Nywi3aPaWSAiPyg8Io+LyB4RyYlI0+ReDiHkUkNVo6r6bwUJQb+IPCYi1YUF0iucGPX9MREJysj8tUBE3vnK3FSYn66TkSfplk+KiIrIs6q6S1U/OBXXQ2YOLqDIVPEtEblfROY556pE5EsyMpmIiLTKyCNtERFR1TIZeS34CidE5PXOuepR/yLOuZPTlDsh5OLlEyKyQkQ2O+cqReT6ws911D7zRn0/X0aeJnXJyNz0DTM3lTvnPmtP4pxrc8592Dk3W0aeqv+rqi6digsiMwMXUGSqiIlIj3MuqaqbROQ3R227V0TeXBChh0TkzwUnry+JyF+r6gIREVVtUNW3TlfihJCLiqCqRl75JyOapISIxAvi8D89yzHvU9XVqhoVkb8QkXsL0oVvysjc9TpV9RfavOEsInRR1XeO+nmvjMgUclNxgWRm4AKKTBW/KyJ/oaoDMqJxuueVDc65XSLyURkRobeKyICIdIhIqrDLF2Tk6dVDheOflhEBOyGEeOUnMrJgeuVftYiUycgTpadF5MGzHPMNEfmaiLSJSERG/rBFnHMnROStIvJ/RKRTRp5I/W85+/+lV4nIM6o6KCPz2cecc0cm6ZpICaDOubH3ImQKKfzlXlxElnGCIYQQciHAJ1BkRlDVNxfEnOUi8nciskNEjs5sVoQQQsj44AKKzBRvFZFThX/LROTdjo9DCSGEXCDwFR4hhBBCiEcm9ARKVW9X1X2qelBVPz1ZSRFCyHTAOYwQcr6c9xOogunYfhG5VURaROQ5EXmPc2735KVHCCFTA+cwQshEmEhpjE0ictA5d1hERFW/IyO6lledfMrLK1xtbbEmbP9gErarKsR2cVdXUwFxIID7d/cmIG6qx/39Ptw/k0VLDp/iA7loRRTidBL3D0Um0n1q4vwY28+90E2lsxCnh7FvI2UhiP1BP8S5DB5vUdN3Q4MpiCMRbM9n7mXepJ81l2vvfTKF+WRMfrk8NmDHSlN9JcTpTN7EGWw/i8fb81VFsf8GEti/tVURc3yxvZDpa/VjfMbvMCa241TNODVdd0bf23t3BiYBu78zfd3WOQBxTSwMcU8/fg47O1q7nHNnqxU203iaw2KxmKurqzvbJkLIRcqxY8dedf6ayApgjqDdfYuM4dVTW1srH//4/zod//zJPbA9GAxCnDP/y77nbVsgbmgog/ib926H+H/eeR3EsSj+J9fZGYe4LILbN1x7BcQn9uF/HPNW2MnUy6LHb+K0ie3bVbvAwvjg8S6Ijz2Pfbvq8vkQVzXXQDzQgcdbgiG8N88+dRjiZctjEJeHccExmMK+6B3CvgoEcSjuO9oBcWt7D8R9/cMQ2wXS//7wLRAfP4X37tipTojbunFB1NaG/XHbRuy/R7cfgPh9b1wJ8anW4gJt3mxczAUqsa9yGbPYzGPc3tEHcTiM41T9ZuEfxr6MhDB2Zpg6s/gORHH/jFmM/82/PwLxO29cDPF3Hsb1x7/8w58dk9LE0xxWV1cnn/nMZ6Y8KUJI6fChD33oVeeviWigzvZr7RkrBlW9S1W3qeq2oaHBCZyOEEImlTHnsNHz18DAwFl2J4RcqkxkAdUiWC9oroz8STrgnPuyc26jc25jeXmF3UwIITPFmHPY6PkrFsMnh4SQS5uJvMJ7TkSWqeoiETkpIu8WrHd2Bvm8k2RqlJ7DvJVKpfEH5UHUVkTC+BopPYSvvdS8Bjqwrxfiy9bgK7eseY1WHsMF3lA/bp+3Al97nYl9ADc6H9TcPPss5rZpk11c2rbsWhfb+/dvvwDxh98wG+Ks4muao7tPQFxTg6+FEsOou1kwC1+XPtGCr9g2rGvE85n3RPuNPq3BvOIbSuD2ZAI1VukE5pM1mqWqMrz3uw5hfr4c7t/Xj6+lutu7MZ8Unv+hl/Ep7pZl+Eq8Ooj994O9xf59dncbbCuz49jhva2cUw/x1Y34H3cqgvd+pNJEkcOteO0LZ9dCnDSv7Pymb6oC5RAf3xuHOJfC16dxo5ESG5cunucwcnEQvhNfM6fuXj1DmZALmfNeQDnnsqr6+yLyMxkR9Hy1UOOMEEJKHs5hhJCJMJEnUOKc+4mMFGokhJALDs5hhJDzhaVcCCGEEEI8MqEnUF5JprOy73BRn5G3BjYOdS7zF6AWpKMHdSv9CdRybFyAGqctW5ogPtWG2g2/Of/SRagV8cXG8M8Zk1f3Vtp0WehVt50d1JUcSaJGp68ddS/hMP5peW83/in8w7tRg/Xe6xdCXFWF+R3uwj/rV7/RaJm+9JnNuUHUqx3pwHtZXo56s8Sw8X0y158z3kgL582C2OqMeuLG9iBrdESK7VVEUH/X2Yo2CrVXoObroZdR51QfK55/cNh4TuVMbDRJ8aOoB/vxYfw95823rYO4wXhe/e0XfwXxn/7BtRCHa1DPtuvZoxA3NaD/2bee3AlxMIz5PPAoWlqkU1ajRUhpQc0TmQz4BIoQQgghxCNcQBFCCCGEeIQLKEIIIYQQj0yrBiqfy8tA/1Dx5Kamly3ZtWg++i75g6hT2bcTfTv/x2+iNiSZRG1JPoO6oUAYtSC+mC3NYv1svK433at8LyLlqNF5fh+6tF+5wpr2YecsiqBOZfYs9O6xteXKwnhtN683fW/qs4m5N3XV2L4Yz6684vWFy1FDVFWJ8bxm9C46fgo1WTlnauGZem0235pq1AFZF61EytTCMz5SORxastbo6X7eajRgAfzodJj8Q9Fi/zlz77IZU3vO+JH5TdtlZthtewbL9BzO4fFf+Ph6iE8O4fZVS7Dvf9wdh7ilFzVM0TDmu3zOHIj3HMG+OdCFejtCCLkY4RMoQgghhBCPcAFFCCGEEOIRLqAIIYQQQjwyvRooJ5IeXQrP6GZ8eVzPhcOomwkGUbO0+ZYFEFeEUReUNe35jLakqhL3t+zej1qO1ctrX2XP8RA0MWp8Vs6JiBfueRiPv+Xqudh6Frf7zFK5IWKvHXcIBDHfIVN3cNMy9EEKB8z1Ga+jsGnPGU3T8BDq0zIZFCXlMtieP4C6nmgUfausx9fwEPpA2etJ5rH92rmoqVrXhz5TzbV4PRmj2dJUcbspdSc5H+qxIuZzkDf3LmW0gScG8VqyGWzvX78fh7gign3zyJPYty9ub4d4yzLUf73t+vkQ+831pHPoW/WGqzdCfP+9QgghFx18AkUIIYQQ4hEuoAghhBBCPMIFFCGEEEKIR6ZVA+VTkVCouGarqUOfp7UrmiF2xhgqEkLtxuM/PwHx7R9ZBHEijbqaXN54AQ3Zml2mPp0Vr4yJdR8qxgc6h2DLsgbU8PzV3Y9D/Oe/dz1m5sfcDu77BcRX3YjePNYnKZE0tef859ZcVTdje0f2tEC8fmkDxGpspJ5/GfVjtQ3oazVsNEgJUz8tZerDqfFSqp+N7QWDmEBPP/pqDSax/XQCY5+51w9u3Q3xW29AjVk4h/snEji2kv5i++Gg0Zf58VqcD3M3ZfxEjSdX1tQBLMvhuVMOr63X+jpV4Mc+Y/zR8kPYvv8MvzbMd93y2RCXl6F2kRBCLkb4BIoQQgghxCNcQBFCCCGEeIQLKEIIIYQQj0yrBmp+c6V88Y9vOx3vODEA2493oO4lFkatyLYXUYfjz6MWI5tHDdJAP+p+rHajPFJtMkStSEce/W1WC2q2zsQY9ozSQC1rsF2N56rPYz2xziMHTdN4bR943WKI/UZHo4qaqaQaPZjRzQwPoyYpeOokxPuOHoa4eTlqoGx9t6zV+Rjfp55+7NtUCs+fM3UIU6b2XnMF1ubzh/B8g53Yv4vnV0N8sg37W0x/xfviEIeMz9WPd3RCHC03Yyta3D+VcGZfbCudNp5XDq81Z3yhAgH8vcflsb2mWAXEQ4PY1/EEni9j9GbHBnD/TT4cu1mjJQwavV0fys8IIeSihE+gCCGEEEI8wgUUIYQQQohHuIAihBBCCPHItGqgEmknO04UtSmJXtzeGEP/mI1rsPbcteswzuRQ+5HIGg2U8bOpjOJ6cW4z+k5Zblh17u1jM1orgrl95eljEL/p2iUQO1sn0HjxqM94TjnUAPmiWDcw1YW6lqpK034UdSzO6FyuWdwEsZHhSNacvzyMQ8s6ZGXzqLtJWh8oo9PJ5TGuqMTry6aNBiuH7b+4swO3m/px775uBcT//pM4xNEIjs3O1laI7fVlMkVNl5EwSd7owcrC2Pc+UzcwYfVhps5g1uHxw0bTJEG8F40RvHnRKHqC/fptyzBfNZosW1jRaA9ra6qEEEIudvgEihBCCCHEI1xAEUIIIYR4hAsoQgghhBCPTKsGSlUk6CvqKTIxPL1PUUfz8v5+PD6EWo8KU2OsuSEKcW096mQqIkYrIqiLCYvVPFlfJ1Ok7AyMEiYxav8y1K289OOdEN/6XtTgWJmJ1djY+mg+o5nyK+pgGhpwfyNbkbDp24zRHIX9xnvIZJRxuL8aXY/ZLMkk5pe3fW28kCJGF1QeNrX8bPspk4/pUHtnW5PokzV3Ltbaa6zE66mvwnwypkMTyeL50+bifQFsK5Uzei9To9EXsPouzL7SeFD5A+ZjncW+PH4K9XA3v205xHnT90EfXqvf/N6VNZ0/2IceWYQQcjHCJ1CEEEIIIR7hAooQQgghxCNcQBFCCCGEeGRaNVD5vMhQsqjncEYbEqrH/RtSmF57CrUZA1nUimTaMba6GKlC7Uj8cBvEqX3o7VNZjjXFGhpnQ9y4YBbEflPf7j9/tPv09x94F/o8ZVN9EPuMKOfFbVj3b+16PFcgiLXbjLWQZBO2IBnuoEYFlEyhJslvrX5Mgta76IntqHtprq2EeGAANUZpc75MylyA0cPNa0ZvoWgEx0a8H9tPneErhV5KRuYjv3hiP8Rvvg59r1q7sf2+JA6ucqNryrtRHejDk1kPKiPlE7+tI2juVSaP15ZOYwM5UzfQ3tvDbaj9e0vz1ZhfD449Z/R2Z1R1NB5gmaqxakYSQsiFD59AEUIIIYR4ZMwFlKp+VVU7VHXnqJ/VqurDqnqg8JW/chJCShLOYYSQqWA8T6C+JiK3m599WkS2OueWicjWQkwIIaXI14RzGCFkkhlTA+Wce0xVF5ofv1VEbih8/3UReVREPjVWW4lMXna1Jk/H+3diPbhkDnUqqy5DX6al1ej9U1aOOiC/1WKYmmEp4/UTLMP6Zvb8qSQePzA0BPHeI3shrgqhX86a2UXtyY6n0Xvno++6HOL/fvIwxEdPoB5r6w6MP/reKyFuGcZbuXEenu/5J1ETtfRyrCvojIgqb+qrqdEkGZslWTQHNUp+c3xHD3p6JY0mKZPDe5MxBeRqKtHjS/3YfjKD+w8Pm/bMvS2rxLF06hjq4crKF0L82N4uiAOmPwZN7UD1FzsokMV9bR1AUbx3ts6gaVoiYRz3l69vhPipp0/i8UbwVRXGvuzbtgPiiqV1mK+pM+hCRj+Xwe31Z7hszSyTOYeR8RG+czfEqbtXz1AmhEwd56uBanLOtYqIFL42jrE/IYSUEpzDCCETYspF5Kp6l6puU9VtQ4MDU306QgiZNEbPXwMDnL8IIUXOdwHVrqrNIiKFrx2vtqNz7svOuY3OuY3lFbFX240QQqaTcc1ho+evWIzzFyGkyPn6QN0vIneIyGcLX+8bz0GNteXy++/efDpWuQq2pzOoW/n5SzinPfMI1o9r6+6BOGDqnd1wLXovzarCy01mkxDPR9sn6TfLy8F+3N/Y8UhnAPNXLep0wmHUIJXFsE7fe67HXBPDaIrVk0Ldy/F21BQNxVHz9Lc/3gdxeQ4vTmvx4o5148XcuM7UBTSylmAOdS+pfhRF9Q4PQzw4gBokX9A06Ef9mM9olmoqy/F8xuepbwj7fiiNx2cxlFtvngvxNw/hWAsaPVvW1KurrEIdkdUtHTvRW9wWRL1W0Fyr2SxBP47ThPGN8hlR1MsvoT4rHMbjo6aO4ELFsVKxqBriVAav1WoLfebelxtPsNTwBfGk5rzmMDI+qHkilwLjsTH4tog8JSIrVLVFVe+UkUnnVlU9ICK3FmJCCCk5OIcRQqaC8fwV3nteZdPNk5wLIYRMOpzDCCFTAZ3ICSGEEEI8Mq218ESwIpsKajPCRnfyxqvKTTwf4s9+cxc2nkRvn6dfwPpm5VHUjjQ3YL225mtRFxMyXkR9Q1i/LphHLUl7GnU/O44WdTlr5uFfSTeX47GZIfSkGkzi2rYshBqfGqOhWrkAr33LFddCvO8kaqbad7VDnJhTDfFwH157tBLzuefpAxAvjKLAdqH5o/CVS1FTtXc/5hMIYt8OlBmvo8vw3hw4gfq3VAL7x9absxqmjgO9EDfPQ41YfQ32p89oqDJ+7J+eYSOICxR1Q+XGb0xNDci08eBKp7GtXA6vpTyKH1uXN7XzsCvkZBdqktasQ5+ngKm9lzdlCYMR1DhljEZKTa2/QAA/V4QQcjHCJ1CEEEIIIR7hAooQQgghxCNcQBFCCCGEeGRaNVDOOcmli/qJnKLYwmoxxGf8ZwS1GB/fvADinz6J9cyyDbi9v+sg5hNH3U1jlfEG+im2d+AwegUN96MOR41fTnKUlqV1H9ays+Y/5UHUNC01mp8rjFdPPos6lLYOU5zO9FUkgmvlVevmQByrR51OwJRryynqXDYsQc1QVydqmvZ34PWsrcY6gs0x1PkEDuD2dW9eCXGZ8XlK9eH+4QAO5bpa1M+Z0nny0EstEL/7GtTX9fWhni1odEe2TmJiCPd3o7yS2o0vUtjci8ow9pX4je+TD8+dNrXpwnjrxJfF7R983SKI/+kRrEHZvBDPFw3h+Xym7p/PbP/lCbz3t16OfUkIIRcjfAJFCCGEEOIRLqAIIYQQQjzCBRQhhBBCiEemVQPl9/mlMlbUzrgcCm2OnERvnoZ61NlE6zDdzDzUuax94zyITz2MNcIGB1FHM29pDcT7d6Cu5oarsB7dDZdVQZxKY3vHTL23R544fPr7xCBqZoJBXLuqD9u6eRN66fR1Qih5Y/YzpwFFPu0DqIMZGsbcDiex730d6EsVMPlFwuiLVCbo07RmKWq2olH09Bo21x8oR9+oBddi39YvwvYGjqB+rb4adTnNs2shrrV1DBWvd39LN8RlEaxtl0hhA3ObcOwN9uP11M7CsRpKFK+/dwj37e7DmoqZPOYWMHotU+JRcsaH6fH9qM27bGkTxM/tj0P8UeOx1WE0TiFTOy+VMWMlhZqpqyK4/7rVeH5CyMwSvnO3p/1Zy3B88AkUIYQQQohHuIAihBBCCPEIF1CEEEIIIR6ZVg1U3uUlkSzqPwYSqLuprUZtRjaJOp7DR1AjVb8IDXBqelEX9Pxzj0NcFUWtR7gWa4IdakNfp5XRaogzWdQRdcVR27K4AfO/t73oM6VmrZo0tdmyDoUv3QdQJ7NhA+pKKowv0bDJrbcPY2MdJGUV2BezAqhrOWV0LgMDqA9LBPH4vjZT2894H6nxvZplagPOn4MaJBlAX6VADLdHs9jXAznU4Zzqi0PsE7ye2eWoWVKjO1q3Dn2o1nbjvf7BthcgXmk0UPlR9ekGs3hs5zD2ZciU0RvIYjJ97dgXVbVGo7QTx9JN62Zhg4r7tzXjYAgZDZP1Y0tncexW1KIe7mCLKRRo9HGEkJllMjVNY+mpLiX9FJ9AEUIIIYR4hAsoQgghhBCPTOsrPBWRwKglW9D8KX1X2pRyMW91mirNn4rHcXsmiq/0oiF8VXHrdQshXlyPVgGxvHktlcDXaNEyfLVSaV4J9g7iq4xEcvSrF1yrBh3GYXNtL7bFIe58CtvGF1Ii/QnMZdN6fI2zvAktHxpq8LVN0vypevwUvg5Nm74o95nXQGF8jdSdNuVIsDnxdeMrvof3t5gdTH+FMP8Vi9C2YPkcM5TnY3mUjjY838GDaGOQFLz+x57aB/HOLnP9xnrgyf34CnPT0mK+MWMBUVeBAzufNXfTlDByzThOcw77ftE87ItsDsep74ySSXi6SAD7OmdeL4sp47N8w1KIM/5TuL+Y98WEkIuGS+kV3VjwCRQhhBBCiEe4gCKEEEII8QgXUIQQQgghHplWDZQTkcwoeUR5LZbvKIuZUi29qCuJ1eHfe6f6UIdyzOiGViyfDfGv/cYbIP7Wf/4M4ivXoa3BQC/qZqINqEWZ1Y9ajz/5yYsQu9EXG0JdSz6PGqRr1plSKEZvlUyjLuWN16+D+Fc/2wHxC08ch7jspsV4fj/+qbkK5lMWMaVMKlEHk89iPgOKOp/QEO6fENRwxbN9EM+qxPY6zZ/+pxLY3o6dePxLe/Be+H04lkJBjKNhvB9LmvDeSg7b27n7KMTdpjSMBHD/R/YU7184iNq8ygj21fIG1GuFzXa/4r0wUj0J4a2TtMm9ffcgxEfxdFJ7BZYs2vCrPRBnV5qSRy9uh/h7jx2CuMxhCSVCSGnhtbTLubjQNVG2L7xcD59AEUIIIYR4hAsoQgghhBCPcAFFCCGEEOKRadVAiRNx+aI+wx9BbUi+D0tWDKZRZxIT1Ez5g2gutGIpikN27EFthyU/jMfnk6iz6ew13j9DqOOZfxXqioYHMX/ww8niuarq0NconUO9lT+FXkF+o1l66EnUncxajl5Ap450QBwU7Bu/8RoaHDYeVgnMd3sr6mjyDo+vq8ahdMPm5RBnHbbfP4BleX7wxGGIG4J4vbaMTjSMa/85NXjvWnowTgrey3QK802askEZh5qsZYuw9Izkz72/jsr3oCmDEx+MQ7w7iX2bNj5N/Sm8F/6Q0XOZj/HBkzgOe0y8fi1ey9xjqFnaVYP3tuXQCYj7jH9bxPhEdfVblzJCSCkxEd2S1QxNpp5KZPI1VZOd32j4BIoQQgghxCNcQBFCCCGEeIQLKEIIIYQQj0yrBiqfd5IcLOo5vvaTfthudS2DXVivbMltqL3YUIFajVAGdSh1tRh/5b+egDho6rf19mM+EeOFlDXt33vfLohbj2G+o7PzKXZ1by/qUp58vhWPNfXO7vrdX4d4RRNqnvZs3w9x0Gh88n7jQyWIGm8hZ2rRRYyX0Y8ewPO97U1YH00C2Fc5QU3Tgy+h7ue6dXj8kV70sXrqqTjEsSRqorasa4A4Yy4wM4i+Ujnj45SPmtp95pORzWEHxUJ4Pek0tp8PFs2ZVjXE8NwOzx2tR23f9VtQP/bdex6COD6MY6d9GHOP+HHsROejUVR3Bj20ftWC7VUYzzKf4r1sLsP2NiydB7E6I+IihFw0eNUoedUgTaVm6WxMRHPFJ1CEEEIIIR4ZcwGlqvNU9RFV3aOqu1T1Y4Wf16rqw6p6oPC1Zqy2CCFkOuH8RQiZKsbzBCorIp9wzq0SkS0i8nuqulpEPi0iW51zy0RkayEmhJBSgvMXIWRKGFMD5ZxrFZHWwvcDqrpHROaIyFtF5IbCbl8XkUdF5FPnaivnVOKpon6iLGzqqxntRLRuFsSt21HbcSqHWpLWloMQd7cfgfjD714LcVcrejFljMapPGqKjIWxiFj/gRaIG5qwvc72ok7HnSE6srXlUNMT8uO5VzWiV8+DPzwK8e/dcRnEP34O+8Zn6qmlU7g9a0q75dOYX3wQ87t8A3psVVeYtXgWfZf+++foJXT0yEmI9203GiLTX+UR9MWqW4R1DivnYy27vt1tED++E+9VrKEC4lURU4+uBuvRBQL4UdEwjsV59agbah1lc5VM25uP4zzegXqwn/30ZYiDoWqIm4MYb16J15I1dQSPdKI2b393HOJB428W8OO137pmBcQV5Xi+IeMhZj4mM8pkzl+EEO/MtK/TVNbq86SBUtWFIrJBRJ4RkabC5PTKJNV4jkMJIWRG4fxFCJlMxr2AUtUKEfm+iHzcOdc/1v6jjrtLVbep6raenp7zyZEQQibEZMxfAwMDYx9ACLlkGNcCSlWDMjL5/Jdz7r8LP25X1ebC9mYR6Tjbsc65LzvnNjrnNtbW1p5tF0IImTIma/6KxWJn24UQcokypgZKVVVE7haRPc65vx+16X4RuUNEPlv4et9YbQ0l8/Ls3qI2xmfWb87ogkx5sTPJY/qdbacgfvftWKsuk0bdSkVltTkf1hwzshuJt6FfTta4KaUHUEjkRmm6nKlFV1aBvkrZITz3625E3cnePZjLa29GDdA/34OeVFdcUQdxwNS+GzZ1BkM+491jfKFs7Tz/MPZlKo2aqF2t2FfJIfylP+THex8oM7Xv8sZHKoOaqu6TqKH64VGspZdP4PXUVOAfWTXPaYJ405XoI/X49mMQv7wbn56qyccnmO/6ucXzRWtRFNRUi/q2WABzbRvEtkI+7Bu/qcN3sgPHjjOmXlUhvDc3LsZrDwewvXLjf5bK4eesw3iY+YN4vg1rr5JSYTLnL0LIzDOVmiavjMdI81oReb+I7FDV7YWf/R8ZmXjuUdU7ReS4iLxzSjIkhJDzh/MXIWRKGM9f4T0uZzyPOM3Nk5sOIYRMHpy/CCFTBZ3ICSGEEEI8Mq218CL+tCyrLvoB1dShLsWfRdHRsydQGzJkdDfqw9hvdCiV1dh+Kol+O1ufQZ3LgmasSTanDn2dmhbhXzq/aznqZq67DHVHf/Jvj5/+PmgkRukh9OqpiKD3jjP6rhcO9kLc2xOHuDqGOpmqCuxLvx8TiAVMnb+s6UtTD+2WWzZjvpW4/93/+TTEZadQHxaLog4o5/D8OaPJikbx+uOdqDnKiKllZ4y21HgR5c3+J/Zhrb2v7sE4YO7HrCjG1bPnQnzzNc0QH99VbO+HT+A4GzK17MrKMdnXX4N1AZ3DsVJVhX1THcJ7cSKO+2eM71R/0niQGa2h68d7X4mXLtk89uVwZxziz33pe0IIIRc7fAJFCCGEEOIRLqAIIYQQQjzCBRQhhBBCiEemVQPl8/skVl3UFamidiMbQm3I5qWoxbD15Ex5MnnTa9B/pr0VtRphI0Qa6kddzakg5vPiHvQaymDJL6lsxPprIUXtSThcFI/c958fhW3/79+h7UwoiGvZpClOV6koVHnfh18H8Rf+8VeYnOlb23lBcz5nvIUkj33/o++8APHm16D+a/PcaohbEthZwRD2fdrUHfSF8Wams5hvpfFSSpn+8Q1h+85cvzOXF0lhfkmTjxOjIzK6n44jqJn65pGjuP+ow3/jzitg21f/9jFMxo/nqqxHn6jGEHqG5RKYS4/xH8tlsC+sNnBeLd77HqOPG+zGe99rauX5bd+GUCu4egF+LsiFz52VH4L47v6vzFAmhJQOfAJFCCGEEOIRLqAIIYQQQjzCBRQhhBBCiEemVQOlquIfJVw62RHHZPzoZdRYj9oKv98YCvtQJ9PXh+vBiNHV5LUC4ve+fS3EDz25G+JZlZjP9v3oxWTrs5WHcf8rNhS9gj73pQcwF2O+kzYapJoI6l5SRseyswP1Yovmou+TtV62nlli9WRGE1VehoVTg02ogzkax3yOnsS+aTIeWsfbsJZc0GGGPlMbL2h0OXkjeFOjgSorQ93QGRov43OVDWJ7UdM/mTTqknwBbH94cAj3z+H1dI/y6XrNEvSM+lIaa9e99y1Y2+lnv9oP8ZuuWgDxqrnoN1ZTkYD42SN4LzZsWA5xfQXmmj3cCfGwYnt+UytvQTnGg+bedcXHKmJJpoLwnTh/TWbNMGqeyHiZynFYavAJFCGEEEKIR7iAIoQQQgjxCBdQhBBCCCEemVYNlHNOMtmiliYzgFqJ5DzUZly9cjHEP/vFyxCvW4TeQP3dqC2prEYdkc9crgugbujGLfiudsfuIxC/fX41tme8mfJpFBZt31vU/Ty9Az2lnKLGZ5lpO2i219egBmffz9GXqWEWHp/JoUbImXpo4TLUa+WN71LC6LvmzEH92M+f3APxbVeuhPi4qfe2ZCnWikv1oCbq0RdQh9Ncg/cuZDy8cibfgNG7+XOm+KDRXIXDuN1nfJ4CUexvI/MRTePY6ejvx+NHaba+/pO9eLC5F3Ob0Tep8RAWn6upQD1Z1xCOu2wGr31RXS3EX/wh1ilM5bDv/p+3bIC4el4TxHMqTB3DJOq/dhzEe+dzqKEi08PFrDWZbC4lnc50M9V9WUr3jk+gCCGEEEI8wgUUIYQQQohHuIAihBBCCPHItGqgstmc9HQWtSIZH2qWtB3Tefa5QxA3VqKOxep2Du5sh3jBCvTLOXUSdSq7BvD8b1yL+1dEUSeUTqD2JBhCnUzSoXfQoqXF9i5fhW2Hw7h2DZlreeEw6kx2nsRznzh8AuK3vxk1SGVB1MGEYqgpOkMjlcfz27p+6Tj21WtWom5nV0sHxJVB7Jt4v9HFGF+ltStRd6PGp2rhshqIew50QfxU6wDES+rs0Maxk0yjr5VT3D9gJFQBxfsViJh734n92VhXzPfZRw/AtqjRdw2lsG8H+/Dih41nWNTou9R8jJtXYF++OYV6tmd3noL4T//tl9gepifhKrzWTA7j9ZfNh/jmKxcJIaUMNU8XLmPdO6uRmuz2R8MnUIQQQgghHuECihBCCCHEI1xAEUIIIYR4ZFo1UIFAQBoaitqcelOwLZ/FH7T3o5fQLx9rgfhjH7sG4uVXog4llcH14bx56GW0php1SScPozdRyIfdkxJTT87W5jMF5sKjvIBUMLd0ymhqjO/QyhUYXxFAPVZmXQPEZRFs78TBVogf2hOHOFqN3kK3b5gN8eXXrYN4+/2PQrx8LmqSwnnU2QwOYV+kjUdWWRn2bSiAfRkMYV//8IFdEL//vVdA/K3/+yDEfY31EN/+eqx7eNNS1Ii9fLQP4v2HUNM1lMDry2RMrbw86pJef2ux/tznv/gL2PaaLcsg3voSeoQtWYx1CGNGT+aMHivnUCP13KOoBRxMxiGuqsXPwVWzUM+WSWHfJ41e7kQLagkfefBFiB9/bKcQQshMMJ36Nj6BIoQQQgjxCBdQhBBCCCEe4QKKEEIIIcQj06qB6htIyU9+efR0/KabsNZdIIw6mHmNqNN5768vh7jrGOpW1m+8HOKdB9Ar6eFH0Vdq/SLUwRxqR21HXTXmkzUapz5Tyy+NoYR9xf1TPluLDffN5VB34oymJuMzxkg+U0cwibcyZzRRv/4a1DgFA7g96Ecvop/e8wTE2/airkb8WAsuYOq7/db7N0H8f//5SYhfsxFr4zVHsX8efuoYxAubUKfz+C/wXgbNSM4M4NhInjwO8UEf3uuDp9B369gp9In644+gJqyrG8/33B7U512xqji26xtQI/SaDXjt33joeYg3NaK+LW+kdwPDgxBHg1irbv581Cy9uB/7NmPGVoXxM8sMo94r6EP93ZLZ+LlZNB/jTBL1YWRyKKUaYKS0mMmxcSmPSz6BIoQQQgjxCBdQhBBCCCEe4QKKEEIIIcQj06qBSqczcvRYsQ5XT08j7hBCrUVFELUaEeMdlPehTuXgXtTFRBzqht5y3QqI/abg2kAr6mZsPTTN4nozHzKipzTmXx4r7j8YNz5Cgsd2mVpzCxZVQRw2Gqr+FObujFAmkcS4Gq1/ZDBhat+FsP35jXit7S2ok1m9DH2W9uyL4wkEdTmf/N2rIf7Tu1H3c+NabC/iMP9wNYrGli9fAvH2Hah3yxm92vpVWK9t7hwce2suQ92OlfH86Kn9EL9p9VKInzZ6u7JAMd+P3bEZtqUzOG6zQzgWHt2NHl63VWNf3v8kbn/HzVh7bqgP2z96CjVTS+ZUQ3yizWr/UHuYTWNn5Iz/mctg/mXl+Dkgk8OlpC0h3pjo2JiIjulSHpd8AkUIIYQQ4pExF1CqGlHVZ1X1JVXdpap/Xvh5rao+rKoHCl9rxmqLEEKmE85fhJCpYjxPoFIicpNzbp2IrBeR21V1i4h8WkS2OueWicjWQkwIIaUE5y9CyJQwpgbKOedE5BURRbDwz4nIW0XkhsLPvy4ij4rIp87VVqw8KDdtmXU6fuqFw7B98/o1EO9qQ03SggbU4URMfbicH3UyYT9qN/ymhlgyi7qjaAy3p7Koo6mfhTXKDh3E2nn+IGpPhgaL7fmC5679lirDc/kVNUADJpegYl+kxehUjCbKyL2kOoZ919ePfVEei0B8882oOeqLY/t5Z+sCYl/6fXhv3Ek0Uoo3Yv/sOdxl8sV83vMb6KXk/zYK1mLlqOOpqcMHDMZaSVq78SeV5Zjvwsa5EL/Ujvf6rz72doiTXZ2nv2/rxZqOrW1xiNt7sa2mBsyl3OTyvttQf9XRjR5WPb0YL5lfDXFZGMdOtdE0pbOoaXLm3kaMNjGdNz5TqdLxgZrM+YuQyeLOyg9BfHf/V2YokxEuZR3TRBiXBkpV/aq6XUQ6RORh59wzItLknGsVESl8bTxHE4QQMiNw/iKETAXjWkA553LOufUiMldENqnqmjEOOY2q3qWq21R1W3//wHmmSQgh58dkzV8DA5y/CCFFPP0VnnMuLiOPum8XkXZVbRYRKXzteJVjvuyc2+ic21hZGTvbLoQQMuVMdP6KxTh/EUKKjKmBUtUGEck45+KqWiYit4jI/yci94vIHSLy2cLX+8ZqK5UVOdxe1Hd88I7rYPvhPagV8SVQl/Ojh9CLJxDFCe3ay7FWXtCPvzFGIqjd6GjrNfujliPgx/2TQ+jlFI1i9wX9qF0pixR1RsGw8ZAyeq05fvTqEYf7h4wPlM+H9c72HUMvn4oybH84gyKobBL71m+uPRxC46hdxxMQX3EZ1sK75bULIf7qvdsgvnL9Aojnz8H8hgbw3s9uroZ4cBD7/kO//zWIg+befvI3sXZdKoH3ZjiF15PPmTqHw7h/g/GhChtNW5fRdCWTRV2TKTsoT2w7CXFTHY7jt70W9WZq7r0z+rjDnaihShvfprBi3/gVNU65DLYXDOC1zTHXLlXoS9VyCK99di3qz2aSyZy/CJksZlrzRCaH8RhpNovI11XVLyNPrO5xzj2gqk+JyD2qeqeIHBeRd05hnoQQcj5w/iKETAnj+Su8l0Vkw1l+3i0iN09FUoQQMhlw/iKETBV0IieEEEII8ci01sIT58RlivqMf/8G1kPz5XE9tzaFupdlZbh9zeYmiMsrUOuxYiF6Be1vQZ3owUPYfnMdajvypl5dOoO6o6yptVcZRq8i9Re3Wy+dbAI1P7k85p41TkXO+EAlUpjLwDBub+9FD60eo9+aXYvX6lMcCicG4hCHBPPfu78N40OdEFf48V7t23UM4vLGaojnGk1RZwL1ab9/x/UQ3/XJr0FcF8P2mpZh7budz6Pmq74J+y/ow/6vm40ar94W7M+M8dnKGR8un6+4/ZEnDsK2mlozzuevhfjq114O8QvPYp0qNXX+llfgtezuMX5medwuOBTE+bG9gGBfzJ+LGq2GRqxbuH5hg8mPEEIufvgEihBCCCHEI1xAEUIIIYR4hAsoQgghhBCPTKsGqqmxUv7nx15/On7p8RbY/q0HUeuRqzI1uCLoL7NqDepcwiHjP2N8mrp3Ye262kqsCeb3oa7laAvqWuqqcL0ZDeHxtjf9o7QkAaOB2nsU65UtWIQ6k1QCdStp4ws1EEcfo2GH+6vR9AwM4vad/egbNase+y5q9GRqfKuc8SYSh3qxlPGZ6olj31svo16jnAkZX6pv/HA7xO//jRsg/ta9j0P8jt9Cn5XfeNtVEG+pmQdxbSPq5URQA9Xag/q5ukrsn4oY9l9La7GW35xZ2NZ9j+C437IBq4jsfglrRFYYDzHrEeabXQVxdDgOsdXXmVErYsaKmjqNNTVYR9CiVvVEERQh5BKAT6AIIYQQQjzCBRQhhBBCiEe4gCKEEEII8cj0+kCJX0SKNdbWXbcStq67boXZ3/jXCOp2HnwKvZTSnahTectbUCOlKdTlZPIo1sgF0Mfpg79zDcQn9x7B+BQa6qSMN1MqWdQp5RW3NTaipsinuJY1llIiSTz+pNFAWd+l5QsrIc6hnEtOxrF+WiZv6/qFIC4vxziRRs2T+vD86aytt4b5+4wGKqPW98rofIzGq7p2NsSpDN5bNflvXIE6I2e8lLQc7+WBY9herAm9j/JDqOkaGsL7UVdZ1ES1tuO+YvzOhoeOQhyMYR3CjK2LGMRx29aL505nsS8rjH9aGm+9lIVwe86hJipktjuTT974odnafYQQcjHCJ1CEEEIIIR7hAooQQgghxCNcQBFCCCGEeGSaNVCWsbQSfhPXQXT71XVmuzMx6nC23LIG4m3PYz23dHvcHI9CpGN7UGN1zduvNfsbHyo4P2pqfnrfryC2llJ5h9eSTBsRk7nUyhg2EDOaJWc0VtEK3J4OYG28wWHsO38I75Wa2nw+cy/VJmg8tnLGF8veab/RRGkA82+MoYYsO4TCnpC5/oDxlVKjQzq2Nw5xVx/WvmuuwXpvoTKTsamNl00Xz/fkzhOw7arVqKcKGL3WoKlrGMVLlYzRh53sxLHlTN/3D+K9tBqqcNj6meG0kMvbe4uYoXXGWCCEkIsRPoEihBBCCPEIF1CEEEIIIR7hAooQQgghxCMzoIGyOqVzbbPrO6utsPvb7fbysKbXxisrzXbbHvpMXfPmlWZ7mYnt+c+oOnYav9WZZFDjlPWhDqbcCGF8RscSshoco1upKMfzOYfn608azZMRtuRS2F7Ah7odzWG+eaPT8QeNdxDufkbPB8KYb8DUa/uvB3ZAHK5G76TFC2ZB7Auins1nfKCqjKRp6aolEJ863A2x3/peGS+k3Cjfqp5BrHt467o5EB9uR/3W4ibjgeXDcdQ3aPva3jtT+64MPbSGUDJ1xqdk2Wz8nFg9WyZvPL3Mzcv7MB9CyIVL+E6sUZu6e/UMZVJ68AkUIYQQQohHuIAihBBCCPEIF1CEEEIIIR6ZZg2UE1S72PVbyMTG++ic+ikR6/t05uWNpbGyoK4mG8DjA2e6F5m4qFU52mlrzeGeWaPx0QyeK55A4Uo+gdvr5mCuNpPBYTw+Z8sMGg1PJIT3ImO2G1smSZgG83nU9TjjuxQ09yJgvIiW1VZD3JHBe7v/V4cwASOq2rxuLsSqqOFqNrUCpdcWH8Trz5v6cMbGSmY3oE/U1seKGq38AO481IHXvmYR+of5/efW/jU2VkG8+wj6k6kf70U6Zca98fAaNmMhVmP6wtS28xv9mgkln6EGipCLBWqeXh0+gSKEEEII8QgXUIQQQgghHuECihBCCCHEI9OsgVLBNZtVTyQhcoLeR3qGJsoy0cs5t89UYEzfJyssKnLkpRchzpu2XRaPdTlTzyyBmiIJoybH+irZ2nVWtOP3oc5Fs3g+69uUTmJ+Ab/xsUrj8SnbnqltFzYaK2du3f7uOMRXbpyH+dyP3krBMI6VFYuqMT/Tv62H8PhIOep2hluxQ4NnCH8w4fbuXtx/lBdTcz16Ut30rs0QP/noCxBHTd/4Q3ju1uMnIbal7HLGgytjjJoqynHs5Iwmyh8w+jcjafKf8TExY5m18AghlwB8AkUIIYQQ4hEuoAghhBBCPMIFFCGEEEKIR6ZVA5XNZqW7q6gVqcvGcIdGo3nyYS26M2vLedVajFVLz2I1WmN1l12PFsUj3/rRy7Clahbuu2wx+hItrMbtCxW315djbiFjzNTbh8KVoPXyqUA9l/pRf2brm/mN8CVjRFfOGk8Z3yUxteOsbiafM0Ibc/7m4TrMx2i0KivxeiqCOFb85vr7h/F8sWrUmOV6sD+CZXi8c6ipCg3j9sd3Hzz9/Ya1qN/6xXOnIK6PoR6tIoS5V0TwWvceiUOcydladRifUSHS1FEsC9s6ieZe5bGv8kbPZk5vby25CGF9NEL4BIoQQgghxDPjXkCpql9VX1TVBwpxrao+rKoHCl9rxmqDEEJmAs5fhJDJxssTqI+JyJ5R8adFZKtzbpmIbC3EhBBSinD+IoRMKuPSQKnqXBF5o4j8tYj8YeHHbxWRGwrff11EHhWRT53zZIGA1NVXj27Z7GF9lKImNl5IZ2iUJvpGEsUbA/14vlilrZc2FsXjc+l+2KIONU1rF6Ee7GAX1jfLJFGH0p7Aa22qwfbmVaCXj8/4NHWn0FNLc9he1gpbjKbJyGwkn8b8jORIAlbzlDe+V8biyxnN1ReeeBLidAIPyJSbe+83tQVN/isWYP25th70hUqZsRTVc+uMUln0geprS5z+fvlN1bBt2IwFfwA/hl0DqL8qj2Cdw6E09l0uZ3ycfOZzodgXCSMt3LAS9WVWj3aGpslqrBTPl3X2czmzTNb8RYqMpXmyGqmJtkdIKTLeFcfnReSTgiuWJudcq4hI4Wvj5KZGCCGTwueF8xchZJIZcwGlqm8SkQ7n3PPncwJVvUtVt6nqts7OrvNpghBCzovJnL8GBgYmOTtCyIXMeJ5AXSsib1HVoyLyHRG5SVW/KSLtqtosIlL42nG2g51zX3bObXTObWxoqJ+ktAkhZFxM2vwVi8XOtgsh5BJlTA2Uc+6PROSPRERU9QYR+V/Oufep6t+KyB0i8tnC1/vGd0ovOiVb+85qpmxsvITOqF1nNVRW3IHakliFV98p3P7Ci6nT34fKUHPz/huvhfjBp7Ae2uoY3pqmpahTGTiKT/O6+lGDk+g3mifTlbetXgzxolnY/p4OFMr0dKBGaMjodJJpW8sPz5c5t7WQ5HJ4fH0t5tPyzD6IKypQH7emGf+Iqq4Jjx/qQd3RgRaMd7TEIV47D/+zzOVQ/xYxQ+dfnjgBcShW1KAFjGeVz2iIsknUDDVV4LlyxnMqm8L9c2o609TCC5g6fiFTPK+q3HxOrN4tj5+bxHAK4jM0UtYDbAaZ/PmLjAevGilqpsiFyERU158VkVtV9YCI3FqICSHkQoDzFyFkQnhyInfOPSojf60izrluEbl58lMihJDJh/MXIWQyoRM5IYQQQohHprUW3plY8YTVMJ1RYG2M9qxmaSwN1Ri7B71qObCB7/34kdPfDyYGYVtldR/Et187G+JdB3D7guXLcf9QA8Sn+lGjNGsu6lbue+IkxFtfOArx/IWooVpWVw1xTQx1NJfPx/P3JVET9cIh1Bj5UqjjGTaaKTXGUdevwf747x+gJsuZsXFkAM/3T//5FMQf3bgU4kQDjr1r5lRDnDc+VD5Ty28ogde75fIqiMNtRe8mn6lT6FBCJAGzPWq0dweP4l9/2Vp2gRx+jNNWY2V8otTWrvOhZ5j9HFpbqUAA+956gmm2tHygSOnhVcM0lmaKmigyE/AJFCGEEEKIR7iAIoQQQgjxCBdQhBBCCCEemQEN1Gh9xFi+Tl41SFZT5W37zlbUMK2ZHx6jvXPn29XRevr7oB81P6EAth2Log7l6itNfbI4apyCtegVtCCGxwfKMJfhDtRULTaaq+tmY3sP/aIT4o7BUxD7TddU1ddCfGhXD8QfeccSiIMBzLejD/P99/u3Q7xlw2UQX2Pquy19+0qIs2b7sNUZmQuID+D1VZSj6WsogPkd645DnDqBOqXFy2cV2+4znlmmrmF5BebS0Yvj8HCn0X/ZYWx/DcpjrtVR/JjHjIeWM5onZzRTti5h0EwbGVM7L3NGjUpCJsZka5zoO0UmAz6BIoQQQgjxCBdQhBBCCCEe4QKKEEIIIcQj06qBig9l5AfPFLU1112+ELY3RLvMEVE5N2OJQcbSRKFW4+WtWLB9zW9fZ/a3mifbPmpXksmi95M/hJqfnNGdZB3m4s/htag5t/XesfXN1Ahl7nzXKmxf8darD/N5+xtx/7/7EvosraxATVfTqgqIP/zxyyH+5/ufxoSNGdF8U/vumtVYn621rxLizDxsLm10OP4gaqCyfuzPygjG2RrUPPlNPbeuXtQh3f9sN8RLw2juFAwVj0+b2nXNTag3E3OvfQHcP2zGVTSG+/clcHvQ4XazWS6fjXUZnfkcGMsrSafxXmWy2Nfqw+2q/L2MlDYT9aEaqz36VF0acKYjhBBCCPEIF1CEEEIIIR7hAooQQgghxCPTqoEaSvTJs7vuPx0/9gzqVGZHYxDPrZgF8ZLZKzBe0QhxXZM3H6hsD/rzuDPsa8byocL2HrivBeJ0qig+ec26BbAtl8K2H997HOLjrZjMqcPo43TXr2Nf1NehXswXRG8h3xl6LdQwOaObyWbQd+oPfht9mJyaemmZBMShxmqI33XbjRDPbcB7/U+f/wHEi9agr9SN63EsxDuwv77xyMsQ11WgZuqtr0GfqKEEapZcFoVCxqZLZtXi/ZjfiO0/t+cIxB0/L/bPdeubcFsn9r0EUb92xSbs65tqMJkju05A3ImWW7K9E+9d2HhgxSpQA5XLmM+F8bwSowcLBLG9lKlz6PIXpw+Urz7pyT+IupeLB6/3cqz9vWqqphpqts4PPoEihBBCCPEIF1CEEEIIIR7hAooQQgghxCPTqoFyeSfpoaI2xIVQR9M6hN463fleiJ86uR/i8A7UUDX70YuoPIoaqVuuRm8iI9OR9965wWYs5wa1Hj9/xngdBYrde/XGObApEkOfo/VzsDZdb1cbxJ94LxofDaXTEJ86jvG+flwbzynDvr78MmwvnTZ6MJ8ZGlYW47DvXRA1QScPd+B2IzA71oKarre9Dd+598ZRk7RkLeb7yEuY0NUh7N8nX9wOcdzUEowmsD9+sfcwxId7sT/8EdQpDQ4OQtzXjvkeHy7WxutYgHqvJ1/Gunu//+bFEO96ZifESR/q29KCPlKVlaiR2rwCc7WSJqPAOlP8l8Kxk3PG98m0l8vj9rzx5LpUYb018mrMdG0/y19+4u8h/uTqr0yovUsFPoEihBBCCPEIF1CEEEIIIR7hAooQQgghxCPTq4ESrLulpkZXwvjjVM5DnVCoymg1MhjHswMQ9xufp22HqyC+bQO2v+cXqLmavxG1KYEB1J6E56BX0bqlqD35wJvWn/6+phz1WRmTe3UVtn3zhmaIY7VYq+3UyTjECxfirSzvxtptiQBqpGYtqoH4uw/sg/hkB2qGrrsM9z/Sgn21aSX2rQtg7b+A4vX5/Xj9GkZvolq0TpKTR9ohXhxFfVowhv07mEX9W1kE85kzFzVbGzMo7FkyiD5Te7vQN6oyimNn4Y0NEH/vJ0UN2+Mtxh+sAXP/h0dQ8zQUx755x3ocZ2uX4b3Y0Ya5R01utTEcl1nj0xQwere80TSZWykZU9svGEA9XPIi9YGaakrNG4hcOIxVi28sPvO5PzQ/mZim6lxcTOOYT6AIIYQQQjzCBRQhhBBCiEem9RWeiJPsqD+J9jl8lTFk/nR97dqlEGdNuRG/D4+vKMdXFS6MrxISyTjEjzyLtglBP75m+uUP8bXWqibcvu9Z/HP053bia6Zbrll++vu84GuOkHmFlVG8tspq/NP1ZApfIS1qxFxCIeyLZ3agTcA7bl0G8U9/hq+Nrl+FrwifF3yNk0jiK0BpwfohD3Tgn/UfOYHb1yzFUjZveOMiiCNmLPQP2D+FP/droXQWX9eunzsfdzCOFK3d2J8netA2YngQrzc5jPvHKvGV4zUbr4X4kccePP397Bo8+YCxUOhwmHsuYkoEHeqE+Mf7cdxqGMfOR1bhtXdl8R1cPoO2B43lOJZ85pVcOm363tge5J2azWPZf1ycXGivJuxrHv4p+8XLWK/4zjV2J2qR4LW9C+lzxCdQhBBCCCEe4QKKEEIIIcQjXEARQgghhHhkWjVQKiLBUfKKbAK1FvOrUPejfkzPlzM1JAzpDLa3yvypelcP6lqcKUeiRpdUHURbhXQKdT29bahFyRidUiBU1J5kstj2jna0XFheg5qaTBav1R8wmiC1OhXc/7YbsfRJfwZ1NpctxD+7T5v8Dp5CDVVdDZYj2XILaqps9Y50AkvTdODpJbQPdT3V12OZnR9sexnbT+L1dbdj3199BZ5vYR2OJRdAXU5yEK+v0m90POaT4cvg8XFTEGU4jbqixbOLtgqNtajN23GyC/etxr5dsRzHwve2odYuGsX2ls5fCHFrP1pMPH6qH+IhU5YmYLSFs6J4/lAENVJrZ6NlRCSG+Q8PJYSUPtQ8XTpMRMc0UYsEr+1bpvt8XuATKEIIIYQQj4zrCZSqHhWRARHJiUjWObdRVWtF5LsislBEjorIu5xzva/WBiGEzAScvwghU4GXJ1A3OufWO+c2FuJPi8hW59wyEdlaiAkhpBTh/EUImVQmooF6q4jcUPj+6yLyqIh86lwHVJaF5fZVS07HOozrt+EspuNyqMvxh1D3kzUlJ7Jmf6uhKouhziWRwP1zKdRuJNLY/p4O9Hnq7kEdUzBsS2IUdTLujCo0eKw/gDoS63FlKr+cUQpF8nhtOWPFkzdGSAGfXTtj/Gs3LsHNzoboLRQuN2V1FPVmfb0ogvrb5/ZC/MYA6mx++3bURP38CXwPvqgOdToDQaz9MpBCnVFFCO/Nzn2tEGeTWPrm2UNxiDetRJ+sQaOBihrPsY7WYimcpQtnwbbGYcx9ltEYPX0Az72kHku3xBxqmFbF8HMRCWGuH71iIcQvnMJrjTgc9z97Acd5Jox9eSSO5/ObkkzrGnAslyie5y/LheRXI3IWLcknZiYPMvl49VaaTF3RRDymzqf9sRjr2jxf+4defdN4n0A5EXlIVZ9X1bsKP2tyzrWKiBS+Nr7q0YQQMnNw/iKETDrjfQJ1rXPulKo2isjDqrp3zCMKFCasu0REGhvrx9ibEEImnUmZv+rmXxBP1ggh08S4nkA5504VvnaIyA9EZJOItKtqs4hI4WvHqxz7ZefcRufcxqpqTkCEkOllsuavyobo2XYhhFyijPkESlXLRcTnnBsofH+biPyFiNwvIneIyGcLX+8bq618TmUwXtTOhIK4flO/9Xky3jw+WwsPtwcVYzW6oHKjNQkGULwRN/XUxNRnixkdTV4xXrOoFuJcrpivkWPJ2sZmiANG02TKpUkkiNdibJskqKY+menKkIlzpl6ZrUvYP4gan1g5+hyp0ThlsnhvyiMYX7EAdT8r5lyG+/tQU2WNmHzbDkP8YD96Ef3Ou9Cjy+/HfDNGU7aoEfM53IK6oKsW49PSoTxe79Jm3L714RchvnxeVTGXAOYSCeK1duaw7UofXvt16+ZA3HYC6w7OnYX/sQ+aun0DGTzfkvlYK2/egjKIt1yPfXWoF32k7v3u0xBXmM/tqRT25UwyqfNXV+SC0z2RS5Op1iWdq+2ZZqLX5uV6xvMKr0lEfqAji5OAiHzLOfegqj4nIveo6p0iclxE3nkeuRJCyFTC+YsQMiWMuYByzh0WkXVn+Xm3iNw8FUkRQshkwPmLEDJV0ImcEEIIIcQj01oLT0TFP0ov4fJGN1OB/jbxQdQkqUNhkM9oRerLMe6Oo79NwGiY+vpRK5I3vk/hAGo7rlw+F+LsAGpXbrh6KeYXLJ7PZUzuAaOnMsXkQgHUNCVNLbYgWvFIyvpEGX2YKh4f8mMcrUUdTNyPmqL9XXgvFpl7sXgLXnvXMdTp9KSHII6FUY9mbamOHzkC8cI3rYX4I1nr+YUN5I2m69Rh9N063GXqHprz92fx3paFcYeg4lhVUzcxNreogcrk8F5Hw6hJipVh20eMP1qlqTW36nULMNloHR5/CP/ILJnA3FqOHoV43gL03LK/Vy2pwfY/9ZE3m/2NYC+Lmqk/+ICQEoD6LTIZTKe+aibwkj+fQBFCCCGEeIQLKEIIIYQQj3ABRQghhBDikenVQKlIbpTYxRldz+Ag+jLVVaBWZCCJ243MR1JGs1QTQ3+cVAJ1LT6jA8qZ4wdyqJF67KkTEHd0o9aj/cctEH/o3dWnv0+fsVY1eis/ipoSaXNxxucpk8XtAdMZRhIlQVPLLm1L6XWjxqkyi33XWI99d+RF1EiVHcfra17QAPG/fG4bxL9322aID/Zi31WZe2f1a1bjJc7EPqNxakbNUnMW9XHpPnPvU3h/murRd8pnvI9mVeD23ChPsqF+9EU62tMH8fKGKohvXTMb4t4BvJY9e3shHkygB2QNpiI+4382OGwGwxlj0/qx2Wkid+7t1tOLlAR/s/scRb3IBc1M+j5d6JqnicAnUIQQQgghHuECihBCCCHEI1xAEUIIIYR4ZFo1UH6fSmWseMrZVaj7KQ+hTiWNkidJOfQqShvfJr+prRcyZknpFF6urZ1nvZSCptZdPoB+OrOjKDZ595uXQOx8xfz8phab2jp/eGpxplZdLo/J+Y3mJ694rfkcHp915z5/Ko/X5jPGTGnjezR7fROeL9sF8bETmM+ma+dBfPjESYivaumE+MTmWRAf+N52iL+P8jP51F2bIFZzvQ8+eAjiN2zCe/et57HBlYtQgzW3Cb2QqsJY327hCtQtSaqoKXt+1ynYNM8U1W5tQY+qwX7su1kNuH/A1oQU1Oqls6hBWjEPx2XTLPToOhOroLO/ZxkTMrP/t584Nkb7pBT45OqvzHQKFw13VqK+7O7+r3jafiFzKWui+ASKEEIIIcQjXEARQgghhHiECyhCCCGEEI9MqwaqLOyXtYtqT8fxOHoPJTOoywkHUXczOID11DqMd9GSpajLqa4sxwQUtR+ZNHoBWW8hv+L6sjaM+dx83UpsPoDak8aKYvcebsPco6YtZzRMsxfXQnzqAHr9BEOYa87UujOl9CSj1jfJ6GjsdkPe+E6p0XQ5o8nypVFTteMZzP+THzE+UPvQq2jny6iJ8m3A+m/vj+L5fcZHK2tq0w3m0YspLTg21sw3vldN6M0UCmF+FWWox8v04/mC5cX7Gy7Dex1AGyipjuH22mrUV4VsoUDjCeb8OO7sON53EDUKq9cul3Njf6+yCj2rkcK+f+InO8don8wEn/ncH0KcunuGErkIGUvTNN2ap0tJhzST8AkUIYQQQohHuIAihBBCCPEIF1CEEEIIIR6ZVg2UE59kpKjXiEZR+2GlHuEQ/qAshjqUpvmovagIoRYkZ7yOysz5Bvrw+FgZaj1ODqK/TlcSfafaj6J30KLFqF3JjpLFhE0tN2fXrsb3qfVIHOIAXrrkTB1BfwBNs7J5vLaIqfOXzaFmx2/SSRkPLjU3J2D0YXmHuhirqHrPG1B3c/w4CoGCYbx3665CX6WWNqz/VleDGqS88cnqP9AN8QcjeG92dqD30pDD/ghVoCaqoqoS4q4+1FR1JTG/fL44VrLGYKxvGDt3djmOy8pKU3fP9H1iEPV0ZeVm3BuJ0qOHsG7h6rV4vjPvlsVOE1YDhXUSY76UEEImj7F8pKbTi4n6qiJ8AkUIIYQQ4hEuoAghhBBCPMIFFCGEEEKIR6ZVAyVOBKQmRrPkMyW20kbXEjY6FiPrkbTRBWUENUtDfajNyGRwe6wChUZdx9sgvqwZvYGuXN8I8WAaE4oPFtsP2Np3Zulqrz2sqJNJ5+wBGGaz+APrE5XJYl/WlOO1dvajR1bE+FTlTG09MZqhU614/OzGGojzDq/H58Oh5zexmHs/x9Sis/iNV1HYj3HySqwn9/w2rE+3fBbqjmYbzZPff249XjqBOqDUcHGsPb8Tz7V6GXp8xcpxXPsCpq6hGdd2e7QM+y5hNFdXr54j58ZqoKzvU/bc21N47f6IrZVHSgFqVy5cxtI8kZmBT6AIIYQQQjzCBRQhhBBCiEe4gCKEEEII8cj0aqAkLz4takPKjFSivRP9bbKmPlzOh1qLyihuDwVwPZgzGqjeOHr/hEKm/SzqdA4a76GVjajDiSdN7T6jG8pni9eayVsfJVM7Loi5JDLGa8eEzupQzHY129X4NvUZoydfAHU4OeNL5TeaIlsL75H9XRB/ZOVSiO9/ZBvEd7xrCybchb5KX/4+vuOfvwA1SisWGY2S6Y7IDvQ+6lphfKZmoZ5NY9j/AykcOxFT6y9k6h76zVgMjtIVbVozF7a1dOE4V3Pz6muxrcpqMzay6IEV8KNGqcLcm6hiriLG5EusL5TVQFmNFMb3/XI/bs5anygyE0ynNxCZXMbyfbL3kpqomYFPoAghhBBCPMIFFCGEEEKIR7iAIoQQQgjxyPTWwss7SSVGeTEZe5lwCLUXQ8NJiOtrUctx8iTqXIKKupW00XL4jdlSFneXY22okYpFcH159eVYny0QQO1IuR8b/MetL5/+/u1XLMN9o+UQh4ypVSqJued8RjOl59Z7pVPGd0qtj5OYGH9gVSw+Y1yVTeP+v3ntfIh7ujshvv5y9CI6vLsV4qom9Gla3IxDc++RdoitBqqtHe/dPT0Yv6e8GeJEG46tBYvqIc4b2U8mi4M178f85s+ehQeEi9sTadQ8Rfe3QOw3+rdUBs+VyBp9lRkravR1GaNxSiSxL47uRS2hL4T6t0d2Yh3BWy+vhnj2oiaIH/r5DojLgtY3ilyIjNbhWA0OmVq89jf1bTMDn0ARQgghhHhkXAsoVa1W1XtVda+q7lHVq1W1VlUfVtUDha81Y7dECCHTC+cvQshUMN4nUF8QkQedcytFZJ2I7BGRT4vIVufcMhHZWogJIaTU4PxFCJl0xtRAqWqliFwvIr8lIuKcS4tIWlXfKiI3FHb7uog8KiKfOldbPp9KRaR4yo6+ftje0Y1eQPXVVueC2oxADrUcbUnUftSUobYj6axZEu7//AmsWXbDPPylNFSG9c+sO44/gO1F8sXz/2J/B2zLDuK+v3bNIogzRgfjfKi3MrIZcXnjiRXE7NJGJ3NGVxhNlRrfo5SRtQSMJktMHDA+Unnre+XDBPo7+yBeuqIB4mUrUKPk8qj5uvK1V0OcNb5MAz2ooeoYxuPnBvGjUFtpfKLM7xovHTwGcWcneoYN9BV1R08fxHt/1dJqiHNDqO0LRPFcfnPu+fMwt4jizWk7hp+rvW0YV0ZwLC1diPkk2/BzcL+JfYqap+UL8HPqc/aTMXNM5vxV6ky279M33/mHxeDuCTVFJhl6fJUG43kCtVhEOkXkP1T1RVX9iqqWi0iTc65VRKTwtfFcjRBCyAzA+YsQMiWMZwEVEJErROSLzrkNIjIkHh53q+pdqrpNVbd198TPL0tCCDk/Jm3+GhgYGPsAQsglw3gWUC0i0uKce6YQ3ysjE1K7qjaLiBS+dpztYOfcl51zG51zG+tqqychZUIIGTeTNn/FYrGz7UIIuUQZUwPlnGtT1ROqusI5t09EbhaR3YV/d4jIZwtf7xurLfX5JDjK/6guZ8yIzHKupgy1FJFh1Hrs6UFtR9sg/ob4+EHUudSH0XtpWR1qTyImn1wQNVTZnNERGRXU80Z7Es8W46FW9KxycQjln7ZiLbm23ehTtOIK1ARVJ1GPde1m9CEKG4+qsoitd4bXEjJ9nz7DysdompwRYRnjpJzRcA0OpiDu7Ue9W96UZ2vJ4PUfO4YaqaEcHi/3Po/ny+P2wQ48/zWb0Zeqrx337zh6GOJkFvPJplHD1XsKdUJNFcX/bDdU4b4nW1HLFwlh32XMg45oEG9Oq/HYChoPMJ+5986HPlS9CbyXz+1F/VbQ+DiljV/ai4dwnKcGsG/f9bY1UipM5vw12XitXzaTOpe/2f2hc27/zOf+EGJqcs6Nvfd/+Ym/93T8J++emC/XWGPvXPdvquvuTbTO33SOvfEaaX5URP5LVUMiclhEfltG/ge+R1XvFJHjIvLOqUmREEImBOcvQsikM64FlHNuu4hsPMummyc1G0IImWQ4fxFCpgI6kRNCCCGEeGR6a+G5vKRTRS1JdTVqP8qiUYiHkqizqahETdLmGPrh5LLoC5UyOphtLaj12HUYdUnXL8caX8EYaqRiYeyuTA61Ik8eOA6xCxS1L2V+vNa8ER3lTB2/5uWocYoPxyE+1Ya6msOPo89RbghzTaVw/+Wz0Vdp9XzsW38Gj28xtfkO9KEmqf9EAmIXNfeiC7dXNWPfxtuxL2sXo45nMI7bo6YuYiSM19ccxOOHq/F6fIK6ndQQatDKjK9VOIxjMVKGsfNh/3X3F4VMjVV4L8OmxqNPsa2uNPbdkT2DEF+zBv3JhkydQ78f+yIZwLF15Ah+Lm5Yh/q6F829XFBt8sdLlesvxzqIK5uqhYyNV63GRHQrE23fapwm0pbI2Jofez67/ydXT64GyKsGaaz8xtwu5+7Psft7anVIk6lz8qppmug493pvJzKW+ASKEEIIIcQjXEARQgghhHiECyhCCCGEEI9MrwYqL5Ie5Z/TYTROoqjd8Ntyaz7UepRVGh+pFOpessOokbp6PtbsunoOak3CERR3+ALYPRURXG9mHZ4/1YbpXLG47vT329pQYyMh1PQEcnixGja6FrPWDS4wmp4c9p2inEzy3Xi+4/2YT0cLGjG1vIA6nbmrsO8SphZd3vhA+Yw3kUSNBiyF975hjtEUGc+tqlrcrmYs5EztO3Q+ElFzr8qMpilh+q+rDzVSkRBqrk6cwLGTLsP+7eko6ogiRjTkq8C+2HLVbIgDxkPrA7cshLh5NmqWfD7sq65+WxPS9O16jH1+PH7tZeaDZ+oavvYqU0jReILFyrGvyOQw1f42E/Xf8YJXTZXFq8bK6/nG0jCNtf9kb59svNzrsfb1Om6mexxPVC93LvgEihBCCCHEI1xAEUIIIYR4hAsoQgghhBCPqDPakCk9mWqniBwTkXoR6Rpj95mklPMr5dxESju/Us5N5OLNb4FzrmHs3Uobzl+TRinnV8q5iZR2fqWcm8gUzF/TuoA6fVLVbc65s5VWKAlKOb9Szk2ktPMr5dxEmN+FQqn3A/M7f0o5N5HSzq+UcxOZmvz4Co8QQgghxCNcQBFCCCGEeGSmFlBfnqHzjpdSzq+UcxMp7fxKOTcR5nehUOr9wPzOn1LOTaS08yvl3ESmIL8Z0UARQgghhFzI8BUeIYQQQohHpnUBpaq3q+o+VT2oqp+eznO/Sj5fVdUOVd056me1qvqwqh4ofK2ZwfzmqeojqrpHVXep6sdKJUdVjajqs6r6UiG3Py+V3EyeflV9UVUfKLX8VPWoqu5Q1e2quq2U8lPValW9V1X3Fsbf1aWS20zCOcxTbiU7fxXyKPk5jPPXeec2LfPXtC2gVNUvIv8iIq8XkdUi8h5VndqiOGPzNRG53fzs0yKy1Tm3TES2FuKZIisin3DOrRKRLSLye4U+K4UcUyJyk3NunYisF5HbVXVLieQ2mo+JyJ5Rcanld6Nzbv2oP68tlfy+ICIPOudWisg6GenDUsltRuAc5plSnr9ELow5jPPX+TE985dzblr+icjVIvKzUfEficgfTdf5z5HXQhHZOSreJyLNhe+bRWTfTOc4Krf7ROTWUstRRKIi8oKIbC6l3ERkbuGDcpOIPFBq91dEjopIvfnZjOcnIpUickQKGslSym0m/3EOm3CeJTl/FfIouTmM89d55zVt89d0vsKbIyInRsUthZ+VGk3OuVYRkcLXxhnOR0REVHWhiGwQkWekRHIsPF7eLiIdIvKwc65kcivweRH5pIjkR/2slPJzIvKQqj6vqncVflYK+S0WkU4R+Y/C64OvqGp5ieQ2k3AOO09Kcf4q5FXKc9jnhfPX+TBt89d0LqD0LD/jnwCOA1WtEJHvi8jHnXP9M53PKzjncs659TLym9ImVV0zwymdRlXfJCIdzrnnZzqXc3Ctc+4KGXkl9Huqev1MJ1QgICJXiMgXnXMbRGRIZv5VQSnAOew8KNX5S6R05zDOXxNi2uav6VxAtYjIvFHxXBE5NY3nHy/tqtosIlL42jGTyahqUEYmn/9yzv134ccllaNzLi4ij8qIFqNUcrtWRN6iqkdF5DsicpOqfrOE8hPn3KnC1w4R+YGIbCqR/FpEpKXw27iIyL0yMiGVQm4zCecwj1wI85dISc5hnL/On2mbv6ZzAfWciCxT1UWqGhKRd4vI/dN4/vFyv4jcUfj+Dhl5bz8jqKqKyN0issc59/ejNs14jqraoKrVhe/LROQWEdlbCrmJiDjn/sg5N9c5t1BGxtovnHPvK5X8VLVcVWOvfC8it4nIzlLIzznXJiInVHVF4Uc3i8juUshthuEc5oFSnr9ESnsO4/x1/kzr/DXN4q43iMh+ETkkIn88ned+lXy+LSKtIpKRkVXrnSJSJyPCvQOFr7UzmN91MvKK4GUR2V7494ZSyFFE1orIi4XcdorInxR+PuO5nSXXG6QowiyJ/GTkPf1LhX+7Xvk8lFB+60VkW+H+/lBEakoltxkeS5zDxp9byc5fhfwuiDmM89d55Tct8xedyAkhhBBCPEInckIIIYQQj3ABRQghhBDiES6gCCGEEEI8wgUUIYQQQohHuIAihBBCCPEIF1CEEEIIIR7hAooQQgghxCNcQBFCCCGEeOT/B7uXhNbznNfdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sample random batch\n",
    "for batch in dataloader:\n",
    "    naip = batch['naip']\n",
    "    labels = batch['labels']\n",
    "    labels_one_hot = F.one_hot(labels.long(), n_labels).transpose(-1,1).squeeze(-1).float()\n",
    "    break\n",
    "    \n",
    "# Visualize\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10,5))\n",
    "ax[0].imshow(naip[0].numpy().transpose([1,2,0])[:,:,:3])\n",
    "ax[0].set_title('Image')\n",
    "ax[1].imshow(labels[0,0,:,:], cmap=enviroatlas_simplified_cmap, vmin=0, vmax=n_labels-1, interpolation='none')\n",
    "ax[1].set_title('Labels')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d77ba75-03d0-414a-bb24-452ab9a6a63f",
   "metadata": {},
   "source": [
    "## Train Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc142c96-33c8-4d86-84b4-589412b2bc95",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Train diffusion model\n",
    "diffusion = GaussianDiffusion(T=1000, schedule='linear')\n",
    "net = UNetModel(image_size=patch_size[0], in_channels=n_labels, out_channels=n_labels, \n",
    "                model_channels=128, num_res_blocks=2, channel_mult=(1,2,3,4),\n",
    "                attention_resolutions=[16,8], num_heads=4).to(device)\n",
    "print('Parameters:', sum([p.numel() for p in net.parameters()]))\n",
    "net.train()\n",
    "opt = torch.optim.Adam(net.parameters(), lr=1e-4)\n",
    "\n",
    "steps = 100000\n",
    "update_every = 20\n",
    "vis_every = 10000\n",
    "save_every = 10000\n",
    "losses = []\n",
    "progress_bar = tqdm.tqdm(range(steps))\n",
    "    \n",
    "for i in progress_bar:\n",
    "    for batch in dataloader:\n",
    "        # Unwrap batch\n",
    "        naip = batch['naip']\n",
    "        labels = batch['labels']\n",
    "        labels = F.one_hot(labels.long(), n_labels).transpose(-1,1).squeeze(-1).float()\n",
    "        \n",
    "    # Sample from diffusion\n",
    "    t = np.random.randint(1, diffusion.T+1, labels.shape[0]).astype(int)\n",
    "    xt, epsilon = diffusion.sample(labels, t)\n",
    "    t = torch.from_numpy(t).float().view(labels.shape[0])\n",
    "    # Denoise\n",
    "    epsilon_pred = net(xt.float().to(device), t.float().to(device))\n",
    "\n",
    "    # Compute loss   \n",
    "    loss = F.mse_loss(epsilon_pred, epsilon.float().to(device))\n",
    "    # Update parameters\n",
    "    opt.zero_grad()\n",
    "    loss.backward()\n",
    "    opt.step()\n",
    "\n",
    "    losses.append(loss.item())\n",
    "    if (i+1) % update_every == 0 or i == 0:\n",
    "        progress_bar.set_postfix({'Loss': np.mean(losses)})\n",
    "        losses = []\n",
    "\n",
    "    # Visualize sample\n",
    "    if (i+1) % vis_every == 0:\n",
    "        with torch.no_grad():\n",
    "            net.eval()\n",
    "            sample = diffusion.inverse(net, shape=(n_labels,64,64))\n",
    "            net.train()\n",
    "\n",
    "        fig, ax = plt.subplots(1, 1, figsize=(5,5))\n",
    "        ax.imshow(labels_to_color(sample[0].cpu().numpy(), False, enviroatlas_simplified_labels, enviroatlas_simplified_label_colors))\n",
    "        ax.set_title('Sample')\n",
    "        plt.show()\n",
    "\n",
    "    # Save model parameters\n",
    "    if (i+1) % save_every == 0:\n",
    "        torch.save(net.state_dict(), f'models/unet_{i+1}.pth')\n",
    "\n",
    "torch.save(net.state_dict(), f'models/unet_{i+1}.pth')\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14eea4ac-6ead-459c-a7c3-304ec0c7394a",
   "metadata": {},
   "source": [
    "## Sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "33f9c9e7-76dc-4e11-aa8a-e2a9f32a4529",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAE/CAYAAAAub/QYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeGklEQVR4nO3de7BdVX0H8O+XSx7yElJIiAZM7SAJdQScW8TBKhqxSFWYduiUGqSSNDMd22JJchNsbiohnYY8UKd2nIkJmgpVAaUwFJU0Sltbi17kIZBgqPKIhASwkYASJPz6x9khZ6+97lnrnLPPPnvd+/3MZM7d++7H7zzuL3v/znrQzCAikqpD+h2AiEg3lMREJGlKYiKSNCUxEUmakpiIJE1JTESSpiQmySL5SZLX9TsO6S8lMekIyXeQ/G+SvyD5c5L/RfJ3+h2XjD+H9jsASQ/JowDcBuDPAdwAYCKA3wWwr59xyfikKzHpxJsAwMy+bGb7zexXZnaHmd1P8rdIfpvksySfIXk9yaMP7EjyUZKLSd5P8gWSG0lOI/kNkntJ/hvJY7JtZ5I0kgtIPklyJ8mFowVF8szs6nAPyftInt3j10FqQElMOvFjAPtJbiL5/gNJJ0MAfw/gdQBmAzgBwCed/f8QwDloJMMPAvgGgE8AOBaNz+RfOdu/G8BJAN4HYCnJ97oBkXw9gH8FsBLAFACLAHyN5HGdP01JgZKYtM3MngPwDgAG4PMAniZ5K8lpZvaImW02s31m9jSAawC8yznEP5jZLjP7GYD/BHCXmd1jZvsA3AzgdGf7K83sBTP7EYAvALjIE9ZcALeb2e1m9oqZbQYwAuC8sp631JOSmHTEzLaa2Z+a2QwAb0bjyuvTJKeS/ArJn5F8DsB1aFxhNdvV9POvPMtHONs/0fTzY9m5XG8AcGF2K7mH5B40Eu30dp+bpEVJTLpmZtsAfBGNZPb3aFyhvcXMjkLjColdnuKEpp9PBPCkZ5snAHzJzI5u+ne4ma3q8txSc0pi0jaSs0guJDkjWz4BjVu8/wFwJIDnAezJ6lSLSzjlMMnDSP42gI8C+Kpnm+sAfJDk75EcIDmZ5NkHYpSxS0lMOrEXwNsA3EXyBTSS1wMAFgK4EsBbAfwCjUL710s4378DeATAFgBrzewOdwMzewLA+Wh8QfA0Gldmi6HP+JhHDYoodUVyJoCfAphgZi/3ORypKf0vJSJJUxITkaTpdlJEktbVlRjJc0k+TPIRkkvLCkpEJFbHV2IkB9DofnIOgB0AfgDgIjN7qLzwRERa62YUizMAPGJmPwEAkl9B4yvuUZPYwGtoE17bxRlFuvS6GdMK657csavlNikXXNxWxu5z9ZnuPv9nJrfc/qkXH2v7mDGm/Dwf6w9/imfMrNAXtpsk9nrku4PsQKPt0KgmvBY44cMHl+m8wr6LQnebIN8nzj1G6FPZbfvyskT89bibFF6vTl6Pqp5/Ge9DIHb3M/W36y4uHGLForW55eVr5zqnCL9A7h1NYYu2P8jlcM+6YtG64D7DzvPft2F2btl9NdZsm184hvu6L183t7BNs/2vFF+fuTfk35cJF8GbLbtJYr53pfCxJLkAwAIAOPTILs4mIuLRTWF/B/J92mbA06fNzNab2aCZDQ4c1sXZREQ8uinsH4pGYX8OgJ+hUdj/EzN7cLR9Jh9PO7HVVWXMrU8VSogj6jJVWmq3ClCa0N2jJ5DhdZfnllcsvKblKdztAWDlImefDm63h9fmj3uVe8wODK9xjrk4f8zhtaOOU9lC/sl5S0mWf4IfGVh3t5kNutt1fDtpZi+T/AsA3wIwAODaVglMRKQXuhpj38xuB3B7SbGIiLRN3Y5EJGmVdjsq1MQ6+Wq/za/Ugb59uy0h/WraEVBKGIE6WsxfXbGZRutjxh63bJ3VxFzFyO2V/LpLBj7lrYnpSkxEkqYkJiJJUxITkaQpiYlI0rpqYtG1TiqogX1qUhuWGDV9s0opjndwkGDdvozASvgyxd1lpac/ZqhPr/sFnNtIFwBWDn0qKh5diYlI0pTERCRpSmIikrT+1sR6oaZ1FhnbyuisXklD1RL+PqIa6gZeEPf33XRU15WYiCRNSUxEkqYkJiJJG3s1MZFxrG8DSbrnDZ24xMB0JSYiSVMSE5GkKYmJSNJUExMpQbslnl5NJDOmJqOJfDK6EhORpCmJiUjSlMREJGlKYiKSNBX2RfpgTBXgeyWys7quxEQkaUpiIpI0JTERSZpqYjKm1KUDtJRAjV1FZDxQEhORpCmJiUjS6lUT890Da+IPaYNqYGkos3apKzERSZqSmIgkLZjESF5LcjfJB5rWTSG5meT27PGY3oYpIuIXUxP7IoDPAvinpnVLAWwxs1Ukl2bLS7qORvWvca2MOsnw2oWBoxQ/ZFctWtfBmdo+zZjVyQCPUe9tWX0nzew/APzcWX0+gE3Zz5sAXBB3OhGRcnVaE5tmZjsBIHucWl5IIiLxet7EguQCAAsA4NAje302ERlvOr0S20VyOgBkj7tH29DM1pvZoJkNDhzW4dlEREbR6ZXYrQAuAbAqe7yltIhERlEs2gPFEnH7XwcMr708t3zVomvaPsZ4KuS7+t3AOKaJxZcBfA/AySR3kJyHRvI6h+R2AOdkyyIilQteiZnZRaP8ak7JsYiItE0t9kUkafXqAC7jWie1lVD9Ktz4VVKnKzERSZqSmIgkTUlMRJKmmpjUBp22VuaUr2KaYrnHcDt3+9qaldIBXPpGV2IikjQlMRFJmpKYiCRNNTEJ6mTQu064NbBOzhk6hupfY4+uxEQkaUpiIpI0JTERSZqSmIgkrVaF/aoKyJ0oc8bi1NTluaooLz66EhORpCmJiUjSlMREJGmV1sSmz5iG4bVzm9bkqy0dTdBQkbrUhUQkT1diIpI0JTERSZqSmIgkrdKamD0zGfs2zHp1edL8bbnfL1tzubsL6Ixy57YVGs/ttyRhvg/qOJ6Atxu6EhORpCmJiUjSlMREJGl96Dt58Mb/pY2nBLeeOO+hlr8vowamuppUTvWv0uhKTESSpiQmIklTEhORpCmJiUjSqi/sNzVedQvovlqnW/z3zeDcLrfBbEwhX8X/NBU/L8V3ricDD8R8uKWl0MxVB+hKTESSFkxiJE8g+R2SW0k+SPKybP0UkptJbs8ej+l9uCIieTFXYi8DWGhmswGcCeBjJE8BsBTAFjM7CcCWbFlEpFK02BvPAzuQtwD4bPbvbDPbSXI6gDvN7ORW+86cOdOWLVvedKz8772htFlLmDjvQc8hWh9k34bZueXV2+a3d1KpjfFUu6zzxDq9sH0d7jazQXd9WzUxkjMBnA7gLgDTzGwnAGSPU0uIU0SkLdFJjOQRAL4G4ONm9lwb+y0gOUJyZO/e5zuJUURkVFFJjOQENBLY9Wb29Wz1ruw2Etnjbt++ZrbezAbNbPDII48oI2YRkVcF24mxMSrhRgBbzay5Qc2tAC4BsCp7vCXqjK3KUzE3+YVt8hus9LT5iaq9laxf9YqYEmK7cfjb5rV+YwqDV3oCC70PndS3alMTqqCdWG2fK1Bpu7iYxq5nAbgYwI9I3put+wQayesGkvMAPA7gwp5EKCLSQjCJmdl3MXpenVNuOCIi7VGLfRFJWh8GRWyhg6LHpPnOoImLSoumJbdOdNXCiP6YgTpB4ZhOXSlGzEsY7n/afm9Sd49la/OTvvhqlcW43IliWtfZ6sSNvVDviynwOYXDlTV+vjl97heqKzERSZqSmIgkTUlMRJKmJCYiSau+sN9U0LRAw1WfSfO35pZjBrSronHr0Ckbguec5Mzc5GsAGlJGB+fwMdoPrLhHfo1b6PcJ1cJ9X0jUpdhf+Bw6wQ+vycduDL/q7pcFPRm8cQzQlZiIJE1JTESSpiQmIkmrtCbGY1/EpKZBC69aHNMAsvuJQdoV13k7v+bqh/IDKRYbbnpmPHeKYqu3zmv168ZZnUA6aWe4IqJzdrtCjTs76dvv1pHqrN2XMG77cT67iCYKEZHxQElMRJKmJCYiSWt7opBuTD6eduKHm8/uBOPZx40uVCMro91QJ3EUBgH0HGOZs88+t0bmcNuV+YQGH4x5e4M1sYhB78r4GBVid36/3FMjc7fpRbsx971esbB4jtBrGFPbDb0Nbi1zzHPe3O3XlDBRiIhI3SiJiUjSlMREJGmVthObPmMalq09WBQro32Se+Pc0Vwjge2B9ms+/s3za92al1sj89XMJnkmB25laPaG4Dart3YwWXAPSqmh19j36370nVxySvE1dWN3X1P3vfR99ideGq6BVqFYv8s/ucr6cEbmB12JiUjSlMREJGlKYiKSNCUxEUla5YMislU131fZdTYPFXKXexoVhhpEFkLyhLHEKZAT+cEZQ8XQxnnbK4j6G0jmg+2kceukefnY3Rmiyhh4sRf6NQCie96h2RsL27jvQ+ELlYg3pvi65/fxDSpQjXp3RNeVmIgkTUlMRJKmJCYiSau+A/jcg8tu/SouktBW4fv3UG3FV/MoNlTd1vL3vtd1pTMIZBm1p+Cs4RGtfzupq0lesQaWX7zaafy6fF24Q7j7XsbVxMqoX7VuHl5ZbVIdwEVkPFASE5GkKYmJSNL6Oiji8Lru7/HLqCsVJvCIqivkzzTR6cztqxvUpf1VX+KIGFgxZe7Tc9sVuoXG1dv+rHCM0GS5bu1ymWeQyHbf25gBE0IDfvbs86OamIiMB0piIpK0YBIjOZnk90neR/JBkldm66eQ3Exye/Z4TO/DFRHJi+k7uQ/Ae8zseZITAHyX5DcA/AGALWa2iuRSAEsBLAkeremG2r3n76RvWBn34249q1gjG1v6UotLuP5VqAF5ngudF7XQP9UxjOJnPdSO0H3jVi4u1l2LEw63bvPl9s/0bdO3doRlDYpoDc9nixOyfwbgfACbsvWbAFzQbowiIt2KqomRHCB5L4DdADab2V0AppnZTgDIHqf2LEoRkVFEJTEz229mpwGYAeAMkm+OPQHJBSRHSI7s/2WHUYqIjKKtbyfNbA+AOwGcC2AXyekAkD3uHmWf9WY2aGaDA4d1F6yIiCvY2JXkcQB+bWZ7SL4GwB0ArgbwLgDPNhX2p5jZUKtjhWYA9ykW+9liaTTtNc+LKXbu2zC75Sl8s3dXNktMKtqdhqomfF9AFZ9K6ycT04k6NBuY70/XHVTBnTW8k5fYnbm+X7MfbV/nb+wa8+3kdACbSA6gceV2g5ndRvJ7AG4gOQ/A4wAuLDViEZEIwSRmZvcDON2z/lkAc3oRlIhILLXYF5Gk9XVQxF4MCui/5+9kzu/29ijU0TyHdBszpiKmk3DSelCbC9VyfUIT2Cyele9U7mtQ656lUBOLaLha/l9LZ9zYHlEHcBEZi5TERCRpSmIikrTKJ89tvqH2daQNCbadiVgbasPjqwrsu9bpFO7csBcGRfR0zq0L/6S8TZzntnLIU8sbS0WxEmpghfZZ7vvvvF7B9wC+epXbVrE4UMGk+Q+2POayNa0HXgQ8JULn9Ql1TB9lVdtCf+sH6EpMRJKmJCYiSVMSE5GkVV8T67L+EGrWtm/j7MK6iU57GrdNlxvSS55jFCbHbR1G38S0vQu1R3LFtCUaT5Oe+F6uUHssV8xEMgXuG+E5yYueOlkrvtpcIbZwE8i+0pWYiCRNSUxEkqYkJiJJUxITkaRVX9hv0ouirK846jZuvfjGtbnlDc9tDB+3EKw7MFwHjVsDHY+9g96tyxdi3fNG9WV2Vroz5LhffMQ0iOyXfsQRc87QF1AxnerdovuKhfNzy0OzP184xhpnZvFi4T5/Ft/MXkPO7OWrt84vbFOF2LEpdCUmIklTEhORpCmJiUjS+loT66Se4damYmpRw8jXBTb8In/P7xYo6Ils4qWtZ3TGomAYwZmkh9e0PwN6aNZ0Xz3LPY/7bDWhSWvVDRLpNMp2TuzWv3yKseavW9yG4ECxsXdMZ/WQTmrGse3idSUmIklTEhORpCmJiUjS+loT60ShD6y7gedG2r0fXzI73y7MnSxl4nxf/av72SSCdZOYJxPeKX/OiGJNR23cxrGY+leoDWTMMcqoTb7o1LdYmCmkuI87wKfLfW77PG3NVm+bFxNeaxoUUUTGAyUxEUmakpiIJK3SyXPfOHi8rRi5+NXlFYvWtti6PEPupKNOzct9DQp1g8ZWznLrqkcZ9Qxf+5zigI2tCwe++kZhUETn93XpFyndK9R/nXc7pm/tpEtbTz7SyaQ4UWnH2Wa7Js8VkbFISUxEkqYkJiJJUxITkaRVWtifOXOmDQ8vf3W5cGZvKPmVoQHahtd5OquGWsiGpjyO4A4u53tdV29rHXtEO0QsmRUq1LqzMnni6NMgd9J77pdBoS+CfJ8x99M/yflyqDCzk+cY7X0NNgoV9kVkPIhOYiQHSN5D8rZseQrJzSS3Z4/H9C5MERG/dq7ELgPQ3MBqKYAtZnYSgC3ZsohIpaI6gJOcAeD3AfwdgAMj6p0P4Ozs500A7gSwJHSs5tvc0CCBvq3cSQyKu3g6b0dMwNHsJc8syocc4tQS3INE1NWWOQMYrlyYbxAbVSdwB3B04rh6a77jbQflvTElphYzlhr7ltGZ362rhV6PDtqtxim5A/inAQwBeKVp3TQz2wkA2ePU+OhERMoRTGIkPwBgt5nd3ckJSC4gOUJyZO/e5zs5hIjIqGJuJ88C8CGS5wGYDOAoktcB2EVyupntJDkdwG7fzma2HsB6oNHEoqS4RUQARCQxM7sCwBUAQPJsAIvMbC7JNQAuAbAqe7wldKynXnwMq5tqNm7n1ChOkSc0gBtQnDw31DZu/lHFAd2u3RuIleFBEw9xL3wLE5Tk+aKceGng+boTlui/jaCx/BKFxtn0/Sn0Y5DMbiZf6aad2CoA55DcDuCcbFlEpFJtDU9tZnei8S0kzOxZAHPKD0lEJJ5a7ItI0mo2UUj4ztjtxxVTRyrUwJy62kduDNcAQrW3fU4ftfmv9UyUcGN+8Srn1zFdON2+kG69z23jEzW5cAf71FVd23z1K664/sn9101YuhITkaQpiYlI0pTERCRpSmIikrQ+F/ZbD+gHFAvqHRVInZ0OCe0VUVB3F904v2T5zt4AcPGNzgxIgSezbE3xGCG+QRBDUi7ku0LPvptGld2oaT19TNCVmIgkTUlMRJKmJCYiSau0Jnb85DdgySkHJwpxe5/6GpSuDNVrIjq0FgZ5czb60oV/3fqgQHHExtCAbZ7fv+n7gV2cfcqYRXy8ixlXr4yGqMNr3fpl/qhjqe5YhjIb/+pKTESSpiQmIklTEhORpFXeTsyaRumfND9fAwvWv4DgzbSvBhKqR8R0gC5OSupMlusEEhy80KPCeYzHjZiXtNDRPtSM0PMhc+uX7ufFrZmN93pnmR91XYmJSNKUxEQkaUpiIpI0hibNKNPk42knzm1vn3bbk0TNtxk6qK+ZmOpVSRgKTD4T0z+3k3qVWycbXtN6Alq1G2vf9nW428wG3fW6EhORpCmJiUjSlMREJGlKYiKStMobu+YKoK3HGRx1XStR24caM0YcpK6z6pQi4ouOuhiatSG/wvkGZtL8rflfe96pMhqeul/8qJBfHV2JiUjSlMREJGlKYiKStOo7gPe4eOQt3zgrFzt1lNVb5+eWfSGO6RqYq8Y1sAInVve9xKKuD9nRe60aWHV0JSYiSVMSE5GkKYmJSNL62gE8pvbQdn3CU88ZmtW6U7Br9dZ5bW0v9dHuAIeRm+S4Ax76jjLeBz3sBXUAF5ExSUlMRJIW1cSC5KMA9gLYD+BlMxskOQXAVwHMBPAogD8ys//rTZgiIn7ttBN7t5k907S8FMAWM1tFcmm2vKSdk8fUItqu2Hl2CJZJQh3ffAepgO+UY7l9WsyELSFVlHhjJpKR6nRzO3k+gE3Zz5sAXNB1NCIibYpNYgbgDpJ3k1yQrZtmZjsBIHuc6tuR5AKSIyRH9v+y+4BFRJrF3k6eZWZPkpwKYDPJbbEnMLP1ANYDjSYWHcQoIjKqqCsxM3sye9wN4GYAZwDYRXI6AGSPu3sVpIjIaIKNXUkeDuAQM9ub/bwZwAoAcwA821TYn2JmQ62O1Ulj15CYYyxxZsBxt6lL49Zx1cm8IsvXtp51yEedt+tptMauMbeT0wDczEZT6EMB/LOZfZPkDwDcQHIegMcBXFhmwCIiMYJJzMx+AuBUz/pn0bgaExHpG7XYF5GkVT4oYrOqaj5uzctqOuhfXWpgdanNFRuQFiNxS7orF+c7XhcnBim++aqBpU1XYiKSNCUxEUmakpiIJK36mlhziaKE2lSv2paNZ3V5PVYuzteqYjp3h+povvpXqp+HMT9AQOST0ZWYiCRNSUxEkqYkJiJJq7QmNn3GNAyvO9h5sqr2OcVBECs5bTJCNSHfgH+9eO8KcUS8T+4+blzDay8PHiP4+VC7wv6IfN11JSYiSVMSE5GkKYmJSNKUxEQkaX3tAC714BbuVzjF8b59AVPCPiudmbijGojWtJA/3sR+HnQlJiJJUxITkaQpiYlI0oIThZRp8vG0Ez/cfPbKTt1Sqh2AvRJpqDnejfnO22VwXpDt1/gnCtGVmIgkTUlMRJKmJCYiSau+nZhqNL2l1zcJqn+FqZ2YiIwLSmIikjQlMRFJ2pjrO0lPTchtChcqG6kNj0jFPH9gvr9lH12JiUjSlMREJGlKYiKSNCUxEUlavQr7vup5m403Y/qzq0gvUp5SvgjropG2rsREJGlRSYzk0SRvIrmN5FaSbyc5heRmktuzx2N6HayIiCv2SuwzAL5pZrMAnApgK4ClALaY2UkAtmTLIiKVCiYxkkcBeCeAjQBgZi+Z2R4A5wPYlG22CcAFXUdDzz8RqTXz/KtSzJXYGwE8DeALJO8huYHk4QCmmdlOAMgep/YwThERr5gkdiiAtwL4nJmdDuAFtHHrSHIByRGSI/t/2WGUIiKjiEliOwDsMLO7suWb0Ehqu0hOB4DscbdvZzNbb2aDZjY4cFgZIYuIHBRMYmb2FIAnSJ6crZoD4CEAtwK4JFt3CYBbehKhiIxLZvl/o4lt7PqXAK4nORHATwB8FI0EeAPJeQAeB3BhdyGLiLQvKomZ2b0AClMloXFVJiLSN2qxLyJJq77vZLuNSBJtK6aBFccZTVrclW7m8NaVmIgkTUlMRJKmJCYiSVMSE5Gk9XdQxDFc/KxzEb8w+1PgfShhrMqkRdXsnZWq87cnZpay0ehKTESSpiQmIklTEhORpNG6aWXW7snIpwE8BuBYAM9UduLupBKr4ixfKrGmEifQXaxvMLPj3JWVJrFXT0qOmJmvL2btpBKr4ixfKrGmEifQm1h1OykiSVMSE5Gk9SuJre/TeTuRSqyKs3ypxJpKnEAPYu1LTUxEpCy6nRSRpFWexEieS/Jhko+QrM2EuySvJbmb5ANN62o3yznJE0h+J5uJ/UGSl9U41skkv0/yvizWK+saKwCQHMimJbwtW65rnI+S/BHJe0mOZOtqFyvJo0neRHJb9nl9ey/irDSJkRwA8I8A3g/gFAAXkTylyhha+CKAc511dZzl/GUAC81sNoAzAXwsew3rGOs+AO8xs1MBnAbgXJJnop6xAsBlaMxuf0Bd4wSAd5vZaU3NFeoY62cAfNPMZgE4FY3Xtvw4zayyfwDeDuBbTctXALiiyhgC8c0E8EDT8sMApmc/TwfwcL9j9MR8C4Bz6h4rgMMA/BDA2+oYK4AZ2R/VewDcVuf3H8CjAI511tUqVgBHAfgpsrp7L+Os+nby9QCeaFreka2rq1rPck5yJoDTAdyFmsaa3aLdi8a8pJutMX9pHWP9NIAhAK80ratjnEBjkIw7SN5NckG2rm6xvhHA0wC+kN2ibyB5OHoQZ9VJTEPPl4TkEQC+BuDjZvZcv+MZjZntN7PT0LjSOYPkm/scUgHJDwDYbWZ39zuWSGeZ2VvRKMt8jOQ7+x2Qx6FoTLL9OTM7HcAL6NEtbtVJbAeAE5qWZwB4suIY2hE1y3nVSE5AI4Fdb2Zfz1bXMtYDzGwPgDvRqDvWLdazAHyI5KMAvgLgPSSvQ/3iBACY2ZPZ424ANwM4A/WLdQeAHdmVNwDchEZSKz3OqpPYDwCcRPI3s4l4/xiNmcTrqnaznJMkgI0AtprZNU2/qmOsx5E8Ovv5NQDeC2AbaharmV1hZjPMbCYan8lvm9lc1CxOACB5OMkjD/wM4H0AHkDNYjWzpwA8QfLkbNUcAA+hF3H2oeB3HoAfA/hfAH/Tz+KjE9eXAewE8Gs0/heZB+A30Cj2bs8ep9QgznegcQt+P4B7s3/n1TTWtwC4J4v1AQDLs/W1i7Up5rNxsLBfuzjRqDXdl/178MDfUE1jPQ3ASPb+/wuAY3oRp1rsi0jS1GJfRJKmJCYiSVMSE5GkKYmJSNKUxEQkaUpiIpI0JTERSZqSmIgk7f8BZ8Ido5m/J2IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "diffusion = GaussianDiffusion(T=1000, schedule='linear')\n",
    "\n",
    "with torch.no_grad():\n",
    "    net.eval()\n",
    "    sample = diffusion.inverse(net, shape=(n_labels,patch_size[0]//4,patch_size[1]//4), device=device)\n",
    "    net.train()\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(5,5))\n",
    "ax.imshow(labels_to_color(sample[0].cpu().numpy(), False, enviroatlas_simplified_labels, enviroatlas_simplified_label_colors))\n",
    "ax.set_title('Sample')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf959ce1-7a64-45ef-8d6a-01dc99b2a57a",
   "metadata": {},
   "source": [
    "## Save/Load Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2429dbd6-3e69-48b5-a1c1-2e00b9f90344",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<All keys matched successfully>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Save/Load model\n",
    "#torch.save(net.state_dict(), f'models/unet_final.pth')\n",
    " \n",
    "net = UNetModel(image_size=patch_size[0], in_channels=n_labels, out_channels=n_labels, \n",
    "                model_channels=128, num_res_blocks=2, channel_mult=(1,2,3,4),\n",
    "                attention_resolutions=[16,8], num_heads=4).to(device)\n",
    "\n",
    "net.train()\n",
    "net.load_state_dict(torch.load('models/unet_100000.pth'))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch_geo_simple",
   "language": "python",
   "name": "conda-env-torch_geo_simple-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
